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SECTION  l.o  GENERAL 


1.1  Purpose  of  t.he  User 1 s  Manual .  The  objective  of  the  user's  manual  (UM)  for  the 
Wind  Factor  Simul at  1  on  MoHel  [WFSM),  USAFETAC  Prolect  1  ? 3  .  Is  to  provide  the 
Information  necpssarv  to  use  the  model  effectively.  This  user's  manual  is 
supplemented  by  IISAFETAC/TN-RO/flDl,  Wind  Factor  Simulation  Model,  Aoril  l°Pn,  which 
provides  an  trr^rfept.h  treatment  of  the  model's  science  and  mathematics.  This  manual 
refers  to  USAFETAC/TN-RD/OOl  freouently.  Readers  may  wish  to  have  the  technical  note 
available  as  suppl  ementarv  reading. 

1.?  Project  References.  The  Wind  Factor  Simulation  Model  (WFSM1  is  embodied  in  a 
constellation  57  computer  subroutines  that,  calculate  simulated  wind  factors  for  a 
laroer  simulation  model  in  which  the  WFSM  resides,  flrioinallv,  the  wfsm  was  written 
to  serve  the  Military  Airlift  Command's  (MAC)  airlift  system  simulation  called 
COLOSSUS.  COLOSSUS  is  MAC'S  attempt  to  test  by  means  of  computer  simulation  its  abil¬ 
ity  to  respond  to  a  conti  naency  at  anv  place  or  any  time.  MAC  perceived  a  need  to 
enhance  the  realism  of  COLOSSUS  hv  addino  simulated  weather.  The  formal  reouest  for 
this  weather  support  is  dated  24  July  1R7R.  This  reouest  called  for  three  weather 
modules:  (a)  terminal  weather  for  takeoff /departure ;  (h)  enroute  refuel  inn 

weather;  and  (c)  enroute  wind  factors.  Efforts  to  meet  reouirement  (a)  beoan  as 
early  as  1077  .  This  UM  documents  software  developed  to  meet  reouirement  (c). 

1.2.1  Project  Reouest.  7WW/D0N  letter  to  USAFETAC/DO,  Reouest  for  Weather  Data, 

? 4  ,ii,i  107R. 

1.2.2  Documentation  on  the  Project.  USAFFTAC/TM-fin/OD]  ,  Wind  Factor  Simulation 

Model T  Model  Descri  ption ,  Anri  1  T?8n.  „ 

1.2.3  Documentation  Concernino  Related  Projects.  Mone.  )  \  ^ 

_____________________  fix-' 

1.7.4  Standards  or  Reference  Documentation.  * 

1.7.4  1  Documentation  Standards  and  Specifications.  Documentation  is  in  accordance 
with  ODD  StaTTdTrd  7435.1-S,  Automated  Data  Systems  Documentation  Standards, 

13  September  I«77,  and  AFM  171-1D0,  300-12,  and  300-15. 

1  .7.4.7  Proorammino  Conventions.  American  National  Standard  X3.°-lRfifi  FORTRAN 
hrnorammino  convention  has  been  adhered  to  except  in  the  followino  cases: 

a.  Extensive  use  of  the  Honeywell  fiDDO-Series  FORTRAN  Execution  Error  Monitor 
( F  x  EM )  has  been  made  to  assist  in  debiiQoino. 

h.  Continuation  statements  are  handled  with  an  ampersand  in  column  1  rather  than 
a  nonblank  character  in  column  fi.  This  permits  runnino  the  |JFSM  in  Honevwell  TSS 
/FORTRAN  as  well  as  in  CARDIN. 

1.7.4. o  nop  or  Federal  Standards.  Documentation  provided  in  the  UM  is  in  accordance 
with  non  standard  74?  r  .  1  -  $ ,  T3  UTeot.ember  1  Q  7  7 . 

1.3  Terms  and  Abbreviations: 

Model:  A  r  enr  e  sen  t.  a  t  i  on  ,  description,  or  imitation  of  a  svstem  or  physical  process 

(such  as  the  atmosphere)  in  another  medium  (such  as  a  computer).  A  model  is  a  sim¬ 
plified,  oeneralited  conceptualization  of  comnlex  reality,  usually  based  on  a  set  of 
simolifvino  assumptions  needed  to  obtain  tractable  solutions.  The  model  is  so  con¬ 
structed  as  to  behave  similarly  to  its  prototype  svstem  or  physical  process  in  some 
sense  considered  critical  to  the  problem  at  hand.  The  model  abstracts  or  preserves 
suitably  chosen  "essential"  nropprties  of  the  system  or  process  heino  modeled. 

Simulation:  A  numerical  technique  by  which  systems  or  processes  are  modeled  dioi- 
tai!y  in  order  to  study  the  behavior  of  the  process  or  system  heino  simulated, 
usually  as  a  function  of  time.  Simulation  has  been  described  as  att.ainino  the 
essence  without  the  reality.  Its  purpose  is  to  permit  drawinn  conclusions  about  the 
real-world  system  or  process  throuqh  use  of  the  simulator  as  a  tool  for  study.  All 
simulations  are  models,  hut  not  all  models  are  simulators. 
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Environmental  Simulation:  A  selectively  realistic  synthesis  of  aerosoace  behavior 
consistent  in  space  and  time.  A  techni oue,  often  involvino  mathema t leal  and 
nrohahil i Stic  models,  of  descrihlnq  or  analvzina  the  environment  of  the  effects  of 
the  environment  on  the  system.  An  environmental  simulation  model  can  stand  alone  or 
can  operate  within  a  laroer  simulation  model  such  as  COLOSSUS. 

COLOSSUS:  Name  given  to  MAC'S  attempt  to  simulate  their  airlift  forces'  capability 
to  respond  to  a  continqencv  anywhere  at  any  time. 

Wind  Factor:  The  difference  between  the  qround  speed  of  an  aircraft  and  its  true 
airspeed. 

Ninety  Percent  Worst:  A  wind  factor  whose  value  is  exceeded  00  percent  of  the  time, 
i.e.,  there  is  a  10-percent  risk  that  one  will  encounter  a  wind  factor  more  neoative 
than  the  00-percent  worst  wind  factor. 

Great  Circle  Oistance  (GCn):  The  distance  between  two  points  computed  alonn  a  oreat 
circle  of  the  earth.  A  oreat  circle  of  the  earth  is  any  circle  around  the  surface 
of  the  earth  whose  center  coincides  with  the  center  of  the  earth.  All  lonoitude 
lines  are  qreat  circles.  The  only  latitude  line  that  is  a  oreat  circle  is  the 
enua  tor . 

Newton's  Iterative  Method:  A  mathematical  technioue  applied  in  subroutine  GRTCIR  to 
solve  a  transcendental  equation  (see  USAFETAC/TN-80/001 1 . 

Sawyer's  Equivalent  Headwind  Technique:  The  mathematical  technioue  used  in  the  ufsm 
to  compute  wind  factors.  For  details  see  USAFETAC/TN-80/r>0l  and  AWS-TR-77-P67,  Guide 
for  Applied  Climatoloov,  November  1077. 

1 . A  Security  and  Privacy.  This  DM  and  the  Wind  Factor  Simulation  Model  (wfsmi  for 
whfeb  it.  is  written  are  unclassified  and  can  be  released  to  the  public.  No  privacy 
restrictions  are  associated  with  the  use  of  this  system  to  include  innut,  output, 
data  base  or  programs. 
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SECTIflM  2.0  SYSTEM  SUMMARY 

2.1  System  Application.  The  wind  Factor  Simulation  Model  (WFSM)  was  desianed  and 
developed  to  he  incorporated  in  the  Military  Airlift  Command’s  simulation  of  airlift 
forces.  This  advanced  computer  simulation  effort  is  called  COLOSSUS. 

Desiqners  of  the  COLOSSUS  simulation  determined  that  three  weather  elements  most 
seriously  impact,  airlift,  svstem  operations:  terminal  weather  for  depart.ure/re- 
covery,  enroute  visibility  for  infliaht  refuel ina,  and  enroute  winds.  The  WFSM  was 
developed  in  response  to  the  third  COLOSSUS  renuirement.  Wind  factors  nencrated  hy 
the  WFSM  enhance  the  effectiveness  of  the  COLOSSUS  simulation  by  enahlino  calcula¬ 
tion  of  realistic  fliaht  times  based  on  climatoloniral  winds. 

The  WFSM  provides  to  COLOSSUS  the  around  speed  of  an  aircraft  on  a  simulated  fliaht 
from  point  "A"  (takeoff)  to  ooint  “R”  (landina),  two  points  anywhere  on  the  nlobe. 
It  does  this  by  calculating  an  overall  enroute  mean  wind  factor  between  these  points 
hy  the  Sawyer  enuivalent  headwind  technioue  and  then  anplvina  the  result  to  a  simple 
alaebraic  expression  for  qround  speed.  For  details  on  how  this  is  accomplished,  see 
USAFETAC/TN-ftO/OOl  . 


There  are  at.  least,  two  other 
appropriate  modifications,  as 
navi  oat  ion  information.  It 
storaoe.  Its  efficiency  is 
operational  wind  factor  model 
required  chances,  can  be  used 
in  number  and  complexity.  It 
weather  information  should  be 
a  part. ial  fulfilment  of  that 


applications  for  the  WFSM.  First,  it  can  he  used,  with 
an  efficient  method  of  aeneratinq  wind  and  areat  circle 
is  efficient,  in  terms  of  computer  run  time  and  core 
especially  evident  in  comparisons  between  the  current 
s  and  t.he  WFSM.  In  addition,  this  module,  aoain  with 
in  other  simulators.  Simulation  efforts  are  increasinn 
follows  that,  a  requirement,  for  a  module  incorporation 
included  in  these  simulation  efforts.  The  WFSM  can  he 
requi  rement. . 


System  Operation.  Before  the  WfSM  is  first  invoked,  a  user-supplied  initial¬ 
izing  roodul  e  must  read  the  wind  data  base  from  a  sequential  file  (card,  tape,  or 
disi';  hut.  as  implemented  on  MAC’S  Honeywell  Series-BOOn  computer,  a  disk  file)  into 
the  computer’s  memory. 

After  the  data  hase  has  been  loaded  into  the  computer's  memory,  the  user  at  any  time 
mav  invoke  the  WFSM  hy  callinq  its  main  subroutine,  EMRWMD.  In  that  call,  the  user 
tells  fnrwmD  the  location  of  point  "A"  (takeoff)  and  point.  "R"  (landinn),  t.he  air¬ 
craft.  altitude  and  airspeed,  and  the  date/time  for  which  the  wind  factor  is 
requested.  The  WFSM  then  responds  by  steopina  the  aircraft,  alonn  a  oreat  circle 
route  between  "A"  and  "B",  calculates  a  route-mean  wind  factor,  adds  that  wind  factor 
to  the  oiven  airspeed,  and  returns  route-mean  oround  speed  to  the  user.  The  process 
is  illustrated  in  Fiqure  1. 


Altitude 

Date/Time 

Airspeed 


User  Call 


FIGURE  1.  Svstem  Operation. 


?.3  System  Conf  iquration.  The  WFSM  was  developed  for  use  on  the  MAC  Honeywell 
SWA-Series  computers  at~Scott  AFP,  Illinois.  In  oarticular,  the  WFSM  executes  on  a 
Honeywell  60S0  computer.  Purina  program  design  and  coding,  the  aoal  of  intercomputer 
compatibility  was  for  t.he  most  part,  adhered  to.  With  exceptions  noted  in  naraoranh 
1. P.4.2,  the  WFSM  should  compile  and  execute  on  any  computer  with  a  FORTRAN  compiler. 

The  wind  data  base  for  the  WFSM  resides  in  the  computer's  core  storage  durino 
execution  and  on  a  seouential  file  (card,  tape,  or  disk)  before  execution.  A 
user-sunpl ied  initial izina  module  external  to  the  WFSM  must  read  that  data  base  from 
the  file  into  the  computer's  memory  before  WFSM  is  first  invoked.  As  implemented  on 
the  MAC  Honeywell  fiOPO,  the  wind  data  hase  resides  on  a  permanent  seouential  disk 
file. 

Con f i oura t i on  necessary  for  a  computer  to  use  WFSM  typically  renuires  a  FORTRAN 
compiler,  linkina  loader,  disk  storaqe,  and  FORTRAN  disk  input/output  capability. 

2.4  System  flroani za t i on .  The  WFSM  is  coded  as  a  constellation  of  seven  FODTPAM 
subprograms  with  the  functions  described  in  Table  1.  The  principal  subroutine  is 
FMRWNO,  with  which  the  user's  proaram  communicates  directly  via  arouments  of  the 
FORTRAN  CALL  statement.  Subroutine  ENRWMD,  in  turn,  calls  other  suhproorams,  some  of 
which  call  still  other  suhnronrams.  The  hierarchy  of  subroutine  calls  is  shown  in 
F  i  a  u  r  e  ? . 


Subprogram 

ENRWMD 


DISTAM 


SPHGLO 


F.PTC  IR 


HOC 


R  P  L  NO 


PPLAT 


TAPLE  1.  WFSM  SUBPROGRAMS. 

Purpose 


Main  enroute  wind  factor  subprogram, 
called  directly  by  user  to  compute  ground 
sneed  from  given  airspeed  in  any  of  three 
modes  for  two  fliqht  levels  and  four 
seasons  of  the  year. 

Called  hv  ENRWMD  and  HDG  or  directly  by 
user  to  compute  oreat  circle  distance 
hetween  anv  two  noint.s  on  tbe  nlnbe. 

Conversion  of  spherical  coordinates  to 
olobal  latitude  and  lonqitude  or  vice- 
versa.  Called  by  GRTCIR. 

Solution  of  oreat  circle  equation  for 
latitude  given  longitude,  or  for  lonoi- 
f.ude  given  latitude.  Called  hv  EMPWND. 

Calculation  of  initial  heading  along  a 
great  circle  route  flown  from  a  given 
origin  to  a  oiven  destination.  Called 
hv  ENRWND  and  GRTCIR. 

For  a  spherical  nrid  system  whose  longi¬ 
tude  boundaries  are  spaced  at  3D°  inter¬ 
vals,  finds  the  lonqitude  grid  values 
bracketing  a  given  arbitrary  lonqitude. 
Called  by  ENRWND. 

For  a  snherical  grid  system  whosg  lati¬ 
tude  boundaries  are  spaced  at  IS  inter¬ 
vals,  finds  the  latitude  grid  values 
hractet.inq  a  aiven  arbitrary  latitude. 

Cal  1 ed  by  ENRWND. 


Retailed  explanation  of  tbe  processing  done  by  each  subprogram  is  available  in 
"SAFFTAC/TM-nn/nni . 

As  shown  in  Figure  3.  thP  WFSM  is  designed  to  o ..erate  as  a  miniature,  subservient 
simulation  within  tbe  user's  overall  simulation  model .  As  such,  WFSM  is  supplied  as 
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FIGURE  ?.  Hierarchy  of  Suh routine  Calls. 


FIGURE  3.  WFSM  as  Server-Model  within  User's  Simulation  Model  . 


a  constellation  of  FORTRAN  subroutines.  Even  the  main  WFSM  module,  called  EWRWND  is 
a  subroutine.  All  commun 1 ca t i on  between  WFSM  and  the  encompa s s  i  no  overall  simulation 
is  done  throuoh  the  arouments  of  the  user's  CALL  statement  and  correspond i no 
arouments  of  the  SURROUTINE  ENRWNO  statement.  0 


?.$  Performance . 


?  .  A  .  1  Capabilities.  The  WFSM  is  capable  of  nroducino  mean  overall  cl  imatol ooical 
wi  nrl  f  ac  tors  for  (Treat  circle  routes.  It  does  so  in  any  of  three  modes:  calm  wind 
case,  ofi-nercent  worst  wind  case,  and  the  mean  wind  case.  The  wind  factor  can  he 
produced  for  two  fliaht  levels  and  for  four  seasons.  In  addition,  the  model,  throuoh 
its  suhroutine  niSTAN,  can  provide  the  oreat  circle  distance  (G(in)  between  any  two 
points  over  the  qlobe.  Since  it  produces  both  wind  factors  and  GCD,  it  is  capable  of 
calculation  both  the  around  speed  and  the  adjusted  flvino  time  between  points  "A"  and 
"n"  if  oiven  the  airspeed. 

?  .  A  .  ?  Assumptions.  The  WFSM  assumes  that  all  routes  are  oreat  circle  routes  or  can 
he"  sTihd  i  v  i  ded  i  n  to  several  lens  each  of  which  is  a  areat  circle  seqment.  It  is 
further  assumed  that  climb  winds  and  descent  winds  pi av  a  neolioible  role  in  deter- 
mininn  the  route-mean  wind  factor.  In  the  model  as  presently  written,  the  simulated 
aircraft  is  always  "at  altitude."  Since  an  arithmetic  rather  than  a  harmonic  mean 
is  used  in  computinn  the  route-mean  wind  factor,  it  is  assumed  that  the  around  speed 
is  less  than  or  equal  to  one  third  of  the  airspeed. 

7.A.3  Limitations.  Presently,  the  WFSM  cannot  produce  a  "simulated"  wind  havina  a 
'TaV-Fto-  day  v  a  r  i  a  b  i  1  i  t  y  .  Hence,  within  anv  aiven  season  and  for  a  particular  alti¬ 
tude,  a  particular  route  will  always  experience  the  same  wind  factor,  renardless  of 
the  nassaoe  of  time.  Requests  for  the  simulated  wind  default  to  the  mean  wina. 

furthermore,  the  present  model  is  incapable  of  producino  a  "forecast"  wind  factor. 
Requests  for  a  forecast  wind  default  to  the  mean  wind. 

Temperature,  aircraft  performance,  fuel  consumption  and  other  elements  considered  by 
tvpical  fliaht  nlannina  models  are  not  included  in  the  present  simulation,  which 
deals  only  with  wind. 

Tnerations  north  of  7F°M  or  south  of  fiO^S  are  not  permitted. 

The  coarseness  of  the  present.  IS0  latitude  by  30°  lonqitude  arid  precludes  use  of  the 
model  for  calculation  of  the  operationally  realistic  wind  factors  unless  channes  are 
made.  Simple  modifications  to  the  exist  i  no  data  base  and  arid  svstem,  alnna  with 
correspond! no  channes  in  the  software,  can  remove  this  limitation. 

Specifically,  the  followinq  restrictions  must  be  adhered  to: 

a.  Ho  route  should  he  flown  directly  over  either  nole. 

b.  Neither  pole  should  serve  as  a  point.  "A"  or  "B." 

c.  P.i  rcumferent  i  al  or  round-robin  flioht.s  in  which  "A"  and  "R“  coincide  will  be 
aborted.  Such  fliohts  should  be  broken  into  smaller  seaments. 

d.  Routes  nr  seoment.s  of  routes  flown  directly  north  nr  directly  south  alnna  a 
lonoitude  line  will  also  abort. 

e.  At  present.,  the  model  allows  flioht.s  onlv  at  altitudes  of  PA.nnn  feet  and 
3 s , n o o  feet. 

f.  All  denart.ure  and  destination  points  must  lie  within  the  latitude  ranoe  from 

'S°*J  to  A  0  d  A  . 

onstrictions  e  and  f  were  imposed  in  order  fo  reduce  the  corp  storane  renuirement.s 
of  t. he  model. 

7  .  c  .  A  Process  inn  Time.  Generation  a  wind  factor  for  a  simulafpd  aircraft,  flvinn  a 
Torn  route  7“i ,  7  ?  A""  nautical  miles)  reouires  approximately  o.i-O.c  .econds  of  central 
procession  unit.  (CPU)  tine  on  a  Honeywell  bORO  neneral  purpose  computer.  The  CPU 

time  depends  on  the  number  of  leas  in  the  simulated  fliqht.  Accord  i  nal  y ,  this  CPU 

time,  roferrino  as  it  does  to  a  very  Iona  route,  represents  almost  a  worst  case  tim- 

i n  n  estimate. 

>  .  A  .  A  Flexibility.  A  number  of  t.he  limitations  of  the  WFSM  can  ho  overcome  hv 
s  eTec  t  Tve  ~mod  i  f  i  c  a  t.  i  on  nf  the  model,  the  data  hase,  or  hnt.h.  Since  he  WFSM  was 


6 


huilt  for  economy  in  run  time  and  storaae,  such  modifications  will  have  the  side  ef¬ 
fect  of  increasinq  run  time  and/or  storaae.  In  particular,  the  resolution  of  the 
WFSM  arid  can  easilv  be  improved  hv  chanqinn  delta-latitudes  and  del ta-1 onoi t.udes  in 
the  model,  alterinq  the  arid  numherinq  scheme,  and  addinq  to  the  data  base.  In  ad¬ 
dition,  the  number  of  vertical  lev el s  and  time  periods  can  be  increased.  With  more 
effort,  WFSM  capabilities  can  be  expanded  to  include  a  "simulated"  (i.e.,  variable) 
wind  by  writino  Part  V  of  subroutine  ENRWND,  which  is  presently  set  to  default  to  the 
mean  wind. 

? .  6 .  6  Error  Detection.  User  inputs  to  thp  WFSM  are  scanned  for  lenal  values.  In 
i terati ve  calculations,  iteration  counters  are  used  to  flaq  computations  that  are 
diveraina  or  fail  i n q  to  converge.  The  Honeywell  fiooo. Series  FORTRAN  Execution  Error 
Monitor  (FXEM)  is  invoked  to  display  to  the  user  any  errors  so  identified.  The  FXFM 
error  nessane,  includino  error  number  and  a  plain  text  diaqnostic,  are  displayed  on 
SYSOUT.  Error  number  bl  has  been  used  for  all  WFSM  errors.  Whenever  FXEM  error  M 
appears,  the  user  knows  that  a  problem  associated  with  the  WFSM  has  heen  detected. 
Each  FXEM  error  has  a  disposition.  Some  permit  the  orooram  to  continue  execution, 
while  others  nroduce  a  nroqram  abort.  Presently,  FXEM  error  61  causes  nroqram  ahort. 
Usprs  can  alter  this  disposition  by  means  of  initializina  call  to  the  FXFM  module. 
Errors  and  problems  detected  and  treated  by  the  WFSM  are  summarized  in  Table  ?. 

? . 6  Data  Base  and  Grid  System 

3.6.1  Data  Rase.  Rasic  data  for  the  WSFM  are  the  IISAFETAC  Sinole  Inteqrated 
Operational  Plan  (SIOP)  winds.  With  a  period  of  record  extendino  from  January  1R7? 
to  December  1076,  the  SIOP  winds  contain  mean  u- component ,  v-component..  and  vector 
standard  deviation,  tabulated  on  a  S-deoree  offset  arid.  This  arid  is  much  finer 
than  that  used  bv  the  WFSM. 

To  prepare  a  data  hase  for  the  WFSM's  16°  latitude  bv  30°  lonoitude  nrid  system, 
simple  averaqinn  is  used.  The  u-component,  the  v-component,  and  the  vector  standard 
deviation  are  averaaed  separately.  The  averaoe  is  performed  hy  axtract.inn  from  the 
SIOP  winds  three  values  1  a  t  i  t.  udi  n  al  1  V  and  six  values  1  onqi  tudinal  1  v.  The  resultinq 
IP  values  are  summed  and  the  sum  divided  bv  IP  to  obtain  the  mean  data  values  for 
each  arid  sector.  In  this  way,  mean  u-component,  v-component,  and  vector  standard 
deviation  are  obtained  for  all  jnp  arid  sectors.  Such  a  data  base  is  constructed  for 
January  (winter),  April  (sprinq),  July  (summer),  and  October  (fall)  for  altitudes 
’6,000  feet  (taken  from  400-mb  winds)  and  36,000  feet  (taken  from  ?SD-mb  winds). 
Usinq  an  in-house  USAFETAC  nroaram  WIND,  the 

A  season  x  ?  altitude  =  0  sets 

of  10fl-sector  wind  data  are  converted  from  u,  v,  and  vector  standard  deviation  to 
direction  (beta  annle,  described  in  Appendix  P,  IISAF ETAC /TM -PO /no  1  )  ,  speed,  and 
variance.  After  conversion,  t.he  winds  are  s  t  o  c  c  H  sequentially  in  a  data  file  in  the 
order  described  in  Table  3.  Each  sector  of  thp  data  file  renui res  one  record  of  disk 
storaqe  cont.aininq  the  three  el  ement s  of  info rmat ion  shown  in  Table  4. 


There  are 

4  seasons  x  ?  altitudes  x  10P  sectors  =  R6A  records 

of  wind  data  in  the  data  file.  To  be  used  by  the  WFSM,  the  data  base  must  be  loaded 
from  the  data  file  into  the  arrays, 

DIR(P70)  Wind  Direction,  Beta  Anqle,  Radians 

SPD(R7D)  Wind  Speed,  Knots 

VAR(R7fl)  Wind  Variance,  Knots'’ 

which  are  located  in  the  COMMON  block  WEA  and  reouire  ? .  6K  words  of  computer  core 
storaqe.  The  COMMON  block  must  be  loaded  with  wind  data  before  the  WFSM  is  first 
executed.  Most  often,  an  initializina  routine  external  to  thp  WFSM  is  used  for  this 
purpose.  That  routine  merely  reads  the  R64  wind  records  in  order  and  stores  them. 


TABLE  2.  ERROR  DETECTION  FEATURES  OF  WFSM. 


Suhrouti ne 

Error  Pescri ption/MESSAGE 

FXEM  Display 

ENRWNP 

Wind  ontion  less  than  0 
or  areater  than  3 

ILLEGAL  WIND  OPTION 

Yes 

ENRWNP 

Latitude  nr  lonoitude  out¬ 
side  leqal  ranqe 

ILLEGAL  LAT/LON 

Yes 

ENPVIMD 

Julian  date  oreater  than 
or  less  than  1 

ILLEGAL  JULIAN  nATE 

Yes 

ENRWNP 

Altitude  index  less  than 

1  or  qreater  than  2 

ILLEGAL  ALTITUDE 

Yes 

ENRWNP 

Subroutine  GRTCIR  twice 
consecutively  failed  to 
converae  when  solvinn  for 
lonnitude  of  Candidate 

Point  #? 

GRTCIR  FAILED  TWICE 

Yes 

ENRWNP 

Number  of  leos  of  the  simu¬ 
lated  fliqht  exceeds  a 
variable  maximum  number 
currently  set  to  SO 

RUNAWAY  PPUTE 

Yes 

ORTCIR 

Subroutine  function  vector 
less  than  1  or  qreater 
than  S 

ILLEGAL  ICON  V 

Yes 

PR  TCI  P 

Direct  north-south  fliqht 
is  attempted 

SOLUTION  NOT  UNIOUE  IN  THETA 

Yes 

GRTC  IP 

Newton's  iterative  method 
fails  to  converoe,  as  indi¬ 
cated  by  an  iteration  count 
exceedinq  a  variable  maximum 
currently  set  at  7.  Because 
of  the  wav  in  which  GPTCIR  is 
used  by  ENRWNO,  failure  to 
converoe  is  sometimes  normal 
behavior.  In  these  cases  no 
error  messaqe  is  displayed 

and  the  nrooram  continues. 

When,  t.he  GRTCIR  function 
flan,  ICONV,  is  set  to  2,  an 
error  messaqe 

SOLUTION  Pin  NOT  CONVERGE 
is  displayed,  and  the  nro- 

No 

qram  aborts. 

Yes 

SPHPLO 

Erroneous  function  f 1  an 

ICONV.  Error  code  IC0MV=10 
returned 

No 

SPHPLO 

Out-nf-bounds  latitude. 

Error  code  IC0NV=1 1 

returned 

No 

D 1 sposi ti on 
Ah  or  t. 


Abort 


Ahort. 


Abort 


Abort 


Ahort 


Abort 


Ahort 


Continue 


Abort 

Continue 


r  o  n  t.  i  n  ii  e 


SPHGLO 

Out-of-bounds  longitude. 

Error  code  ICONV=12 
returned 

No 

Continue 

SPHGLO 

Out-of-hounds  A 

Error  code  IC0MV=13 
returned 

No 

Conti nup 

SPHGLO 

Out-of-bounds  cp 

Error  code  IC0NV=14 
returned 

Mo 

Continue 

Season 

Table  3.  STRUCTURE  OF  WIND  DATA  FILE. 

FI i qht  Level 

Record 

I  Numbers 

Winter 

250 

1  - 

1  no 

350 

l  o  g  - 

216 

Sori nq 

250 

217  - 

324 

350 

325  - 

432 

Summer 

250 

4  33  - 

640 

350 

64 1  - 

6fl  P 

Fall 

250 

640  - 

756 

350 

757  - 

P64 

Table  4.  WIND  DATA  FILE  RECORD  ELEMENTS. 


E 1 ement 

Units 

Col umns 

FORTRAN 

Format 

Uind  Direction 
(Beta  Anole) 

radians 

1-15 

FI  6.  5 

Wind  Sneed 

knots 

16-30 

F 1  5  .  5 

Variance 

knots^ 

31-45 

F 1  5 .  5 

without,  rearrangement,  in  the  COMMON  block  WEA.  Thereafter,  the  wind  data  base  can 
he  accessed  in  terms  of  a  record  number  MRMD  as  follows: 

DIR(NRND)  SPDfNRND)  VAR(HPNO) 

The  record  number  MR NO  is  computed  as  follows: 

MR  MO  =  (Season  Index  -  1)  x  (Total  Number  Altitudes)  x 
(Total  Number  Or i d  Sectors) 

+  (Altitude  Index  -  1)  x  (Total  Number  Grid  Sectors) 

+  (Sector  Number) 

2 . 6 . 2  Grid  System .  The  WFSM  is  required  to  be  economical  in  terms  of  comnuter  core 

storaoe  and  run  tTme;  yet  the  model  must  provide  a  qlobal  wind  factor  capability  to 

the  user.  These  conflicting  requirements  dictated  the  use  of  a  coarse  grid.  The 

original  ’dea,  motivated  by  the  1  a 1 1 tude/1 onoi tude  orientation  of  the  winds  used  as 
input  data,  was  to  have  a  grid  system  in  global  ( 1  at  1 tude/ 1 onoi tude)  coordinates. 
The  qrid  would  have  a  spacing  of  15°  latitude  by  300  longitude  and  would  cover  the 
globe,  requiring  144  orid  sectors  (12  sectors  from  pole  to  pole  times  12  sectors 
around  the  equator).  Closer  investigation  of  the  latitudes  of  the  COLOSSUS  terminals 
and  routes  indicated  that  grid  sectors  above  75°n  and  below  600$  were  unnecessary. 
Since  three  quantities  (wind  direction,  speed,  and  variance)  are  stored  for  each  qrid 
sector  for  four  seasons  and  two  flight  levels,  approximately  0.R7K  words  of  core 

storaqe  could  be  saved  by  reducing  the  number  of  grid  sectors  from  144  to  1  Ofi . 
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The  qrid  system  chosen  for  implementation  Is  in  olohal  ( 1  at i tude/1 onqi tude )  coordi¬ 
nates  with  a  resolution  of  15°  latitude  by  30°  lonaltude.  The  arid  system  extends 
from  7R°N  to  6n°S  latitude  (o  rows)  and  qlohally  in  lonaltude  (1?  columns).  Hence, 
there  are  lOfl  qrid  sectors.  Numherinq  of  the  qrid  sectors  Is  shown  in  Fiqure  4.  It 
is  important  to  note  that  column  1  is  keyed  on  30°W,  not  oo. 


Column  #  — ► 

1  2  3  4  5  5  7  8  9  10  11  12 

Row  #  -I 

1  2 

3 

4 

5 
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8 
9 
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FIGURE  4 .  Grid  System  for  WFSM. 
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A  nrid  svstem  such  as  this,  involvina  undefined  reaions  near  the  poles,  prevents 
fliohts  over  the  ooles  (which  would  he  complicated  by  the  sinnularitv  in  coordinate 
systems  at  the  poles).  The  coarseness  of  the  present  qrid  system  precludes  use  of 
this  model  for  calculation  of  operational  wind  factors.  A  coarse  arid  is  suitable 
for  simulation,  however.  The  qrid  could  be  made  finer  if  improved  meteorol on i cal  res 
o 1 u t i o n  is  desired. 

7.1  General  Descriptions  of  Inputs,  Processing,  Outputs. 

?  .7 .1  Input.  A  user  emplovinq  the  WFSM  to  calculate  a  route-mean  wind  factor  and 
around  snee~d  must  stipulate  the  qlobal  ( 1  a t i t ude/ 1 onn i tud e )  coordinates  of  point  "A" 
(takeoff)  and  point  "R"  (landing),  the  Julian  base  date  of  the  wind  factor  request 
(used  to  determine  season  of  the  year),  aircraft  altitude  (?R,nnn  ft  or  3F,onn  ft), 
wind  option  (calm,  mean,  or  ofl-percent  worst),  and  airspeed. 

To  accommodate  future  qrowth,  the  user  must  also  stipulate  Greenwich  mean  base  time 
of  the  wind  factor  request  and  forecast  hours  ahead  of  the  wind  factor  requested. 
Recause  the  WFSM  does  not  as  yet  have  a  capability  to  aenerate  forecast  wind  factors, 
these  two  inputs  are  ianored  and  may  have  any  value. 

1.1.7  Processina.  The  keystone  subpronram  of  the  WFSM  is  ENRWNO,  which  accomplishes 
most  of  the  aircraft  navioation  except  solution  of  the  equation  of  a  qreat  circle. 
It  also  accomplishes  all  of  the  meteoroloqv  involved  in  calculation  of  the  wind 
factor  (see  Chanter  3,  KSAFETAC/TN-Rn/OOl ) . 

The  ENRWND  subprooram  naviaates  a  simulated  aircraft  over  a  oreat  circle  route  from 
point  "A"  to  point  "R"  in  a  alobal  ( 1  a 1 1 1 ude/l onq i t ude-or i en ted )  grid  system  whose 
resolution  is  IE0  latitude  bv  3Dn  longitude.  The  subroutine  determines  which  grid 
sector  the  aircraft  is  in  and  the  length  of  the  aircraft's  path  through  the  sector. 
Then  the  routine  consults  the  wind  data  base  for  the  particular  nrid  sector,  flight 
level,  and  season  to  obtain  sector-mean  wind  direction,  speed,  and  variance. 
Cross-track  and  along-t.rack  components  of  the  wind  are  computed.  These  are  weighted 
bv  the  length  of  the  aircraft's  track  through  the  sector  and  are  accumulated.  At  the 


end  of  the  simulated  flight,  the  accumulators  are  divided  hv  the  total  qreat  circle 
distance  for  the  fliqht,  producing  a  distance-weighted,  route-mean  cross-track  wind 
comoonent,  an  alnnq-track  component  and  a  variance.  From  these  a  route-mean  wind 
factor  is  computed  usinq  Sawyer's  equivalent  headwind  technique  (see  Chapter  3, 
USAFETAC/TN-fl0/001).  If  the  dO-percent  worst  wind  option  has  been  selected,  that 
wind  factor  is  statistically  adjusted  to  the  10-percent  risk  value.  The  wind  factor 
is  then  added  to  the  airspeed  to  obtain  the  route-mean  airspeed.  If  the  calm  wind 
option  is  requested,  all  these  computations  are  bypassed,  and  a  zero  wind  factor  is 
used.  Presently,  requests  for  the  simulated  wind  option  will  default  to  the  mean 
wind.  In  the  mean  wind  case,  the  wind  factor  is  developed  directly  from  a  mean  wind 
data  base.  In  the  dO-percent  worst  case,  this  mean  wind  factor  is  statistically  ad¬ 
justed  to  a  value  such  that  10  percent  of  the  flights  over  that  route  will  experience 
a  worse  wind  factor  (10-percent  risk). 

2.7.3  Outputs.  The  output  from  the  WFSM  (subroutine  ENRWND)  is  the  variahle  GSPEEn, 
a  real  quantity  representing  the  aircraft’s  ground  speed  in  knots. 


SECTION  3.0  STAFF  FUNCTIONS  RELATED  TO  TECHNICAL  OPERATIONS 

3  1  Initiation  Procedures.  The  Hind  Factor  Simulation  Model  (WFSM)  operates  within 
a  larqer~mode1  .  such  a~s~  MAC 1  s  COLOSSUS.  The  larqer  model  Invokes  the  WFSM  hv  issuinq 
calls  to  its  subroutines,  principally  subroutine  EHRWND.  Refore  the  first  call  to 
WFSM,  the  WFSM  must  have  been  initialized  by  readinq  wind  data  from  the  data  base 
into  the  COMMON  block  WEA.  Users  must  provide  this  special  initializing  routine,  as 
the  WFSM  code  assumes  wind  data  have  already  been  read  into  core.  A  suitable  routine 
is  shown  below: 

COMMON  /WEA/  DIRIA70),  SPD(R70),  VAR ( F> 70 ) 


DO  .11.1  NREC  =  ’  ,R64 

READ  (fc.iii)  DIR(NREC),  SPD(NREC),  VAR ( NREC ) 
iii  FORMAT  (3F15.S) 

.1.1.1  CONTINUE 


In  this  example,  iii  and  i.i.i  are  inteqer  FORTRAN  statement  numbers.  The  file  code 
for  the  wind  data  base  is  f c . 

Job  control  languaae  (JCL)  or  its  equivalent  is  needed  to  eouate  the  loqical  file  fc 
with  a  catalogued  file  name  on  the  particular  computer  system  used.  Assumina  that 
the  wind  data  hase  is  on  a  sequential  disk  file  called  PWIHP.SQ  within  the  suhcata- 
loque  MACRO  in  catalogue  MACRO,  then  the  Honeywell  FOOO-Series  JCL  is 

S  PRMFL  fc ,R, S.MACPO/MACRO/PWIND. SO 

where  fc  is  the  file  code  selected  above.  The  $  PRMFL  JCL  must  appear  behind  the  $ 
EXECUTE  record  in  a  Honevwell  ,ioh.  Structure  and  format  of  the  data  hase  contained 
in  the  file  PWIND.SO  is  outlined  in  Table  3. 

3.?  Staff  Innut  Requirements.  The  WFSM  communicates  with  the  user's  simulation 
model  only  through  the  argument  list  of  the  user's  call  to  subroutine  ENRWND .  * 


Whenever  a  wind  factor  is  desired  in  the  user's  model,  the  user  should  nlace  a  call 
to  ENRWND,  supply  inn  the  values  for  input,  arguments.  Subroutine  ENRWND  then  takes 
control,  calculates  a  route-mean  wind  factor  and  uses  it  to  prepare  the  output  arnu- 
ment,  namelv  the  ground  speed  of  the  aircraft.  The  user's  call  to  subroutine  ENRWND 
should  he  of  the  form, 

CALL  ENRWND  (FRMLAT,  FRMLHfi,  TDLAT ,  TOLNG,  J1ILDAT, 

GMT,  FC  HR  S,  IALT,  IflPTN,  ASPEED,  GSPEED) 

Inputs  to  the  WFSM  are  in  the  form  of  input,  arouments  to  subroutine  ENRWND.  These 
input  arouments  include  the  global  ( 1  a  t  i  tude/1 onqi tude)  coordinates  of  point.  "A" 
( takeoff!  and  point  "R"  (landing).  The  user  specifies  which  one  of  the  four  seasons 
he  needs  hv  stipulating  a  Julian  hase  date  (JULDAT)  from  1  to  3fi*.  The  user  requests 
calm,  mean,  or  RO-percent  worst  winds  by  snecifyinq  a  wind  option  (1DPTN)  from  n  to 
3.  The  user  then  specifies  the  altitude  of  the  aircraft  by  settinq  the  flag  IALT  to 
1  for  2  R , DDP  feet  or  2  for  3  K ,  Clfin  feet,  and  airsneed  in  knots  in  the  input  variable 
ASPEEO . 


*TVe  exception  is  that  user's  may  call  subroutine  DISTAN  dirertlv  if  thP>  need  a 
qreat  circle  distance,  or  function  HDG  if  they  want  an  initial  nreat  circle  heading. 
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To ' accomodate  future  arowth,  the  user  must  also  stipulate  GMT,  the  Greenwich  mean 
base  time  (0.0  -  24.0)  of  the  wind  factor  request,  and  FCHRS ,  a  forecast  hours  ahead 
of  the  wind  factor  reauired.  Presently,  these  inputs  are  ionored  and  may  have  any 
val ue . 

A  description  of  these  input  arauments  is  contained  in  Table  5 . 


TABLE  5.  INPUT  ARGUMENTS  TO  SUBROUTINE  ENRWND. 


Variable  Data 

Name  Tyne 

Definition 

Units/Values 

FRMLAT 

Real 

Latitude  of  point  "A" 
(Takeoff) 

Decimal  degrees 
(75.0°  to  -60.0°) 

FRMLMG 

Real 

Longitude  of  point.  “A" 
(Takeoff) 

Decimal  degrees 
(180.0°  to  -180.0°) 

TOLAT 

Real 

Latitude  of  point  "B" 

( Landi ng) 

Decimal  degrees 
(75.0°  to  -60.0°) 

TO  LUG 

Real 

Longitude  of  point  "R" 
(  L  a  n  d  i  n  a ) 

Decimal  deqrees 
(180.0°  to  -180.0°) 

JULDAT 

Integer 

Julian  base  date  of 
wind  factor  request 

1  -  366 

GMT 

Real 

Greenwich  base  time  of 
wind  factor  request 

0.0  -  24.0  hrs 

FCHRS 

Real 

Forecast  hours  ahead 

>  0.0  hrs 

I  ALT 

Integer 

Altitude  index 

1  =  25,000  ft 
?  =  35,000  ft 

I  OPTN 

Integer 

Wind  option  or  mode 

0  =  Calm  wind 

1  *  Simulated  wind 

(Defaults  to  3) 

2  =  90-percent  worst 

wi  nd 

3  =  Mean  wind 

ASPEED 

Real 

Airspeed 

Knots 

3.2.1 

Input.  Formats 

and  Composition  Rules. 

Input 

arauments  for  subroutine  ENRWND 

are  presenter!  in  Table  5.  Special  format  rules  are 

as  follows: 

a  . 

seconds 

Latitudes  and 
system,  e.g.. 

lonqitudes  must  be  in 
36.5°  rather  than  36° 

dec imal 
30' . 

degrees,  not  the  deorees-minutes- 

h.  South  latitudes  and  east  longitudes  must  be  neoative. 


3.3  Output  Requirements.  Output  argument  for  subroutine  ENRWND  is  the  variable 
GSPEED ,  a  real  quantl ty  representing  the  aircraft's  around  speed  in  knots.  In  addi¬ 
tion,  users  may  call  subroutine  OISTAN  directly  for  a  oreat.  circle  di stance,  or  func¬ 
tion  HOG  for  an  initial  great  circle  heading. 

3.4  Utilization  of  System  Outputs .  The  calllna  program,  e.g.,  a  simulation  such  as 
td A f; ' s  COLOSSUS ,  uses  the  WFSM  to  obtain  the  qround  speed  of  an  aircraft  flying  from 
point  "A"  to  point.  "B”.  In  order  to  provide  the  qround  speed,  the  WFSM  must  calcu¬ 
late  a  wind  factor.  Most  of  the  work  in  this  module  occurs  in  finding  that  wind 
factor.  Once  this  Is  known,  it  is  a  simple  matter  to  calculate  the  ground  speed 
(see  section  ?.l  above  and  USAFETAC/TN-PO/OOl ) .  The  ground  speed  information  is  im¬ 
portant  to  the  user's  simulation  because  from  it  the  adjusted  flyinq  time  can  be 
calculated.  The  user's  simulation  then  has  a  timino  factor  to  use  for  mission  plan¬ 
ning,  war  aaminq,  judging  contingency  response  capabilities,  and  other  applications. 
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3.5  Recovery  and  Error  Correction  Procedures.  To  assist  in  debugging  efforts, 
cons  1 derabl e  use  o7  the  Honey wel 1  Series  6 0 rt 0  FORTRAN  Execution  Error  Monitor  (FXEM) 
has  been  made.  Whenever  FXEM  error  61  appears,  there  is  a  problem  in  the  wind  factor 
routines.  A  plain  text  diagnostic,  Indicating  the  nature  of  the  problem,  will  appear 
on  SYSOUT.  Presently,  FXEM  error  61  Is  set  to  abort  execution  of  the  main  program. 
Bv  a  call  to  one  of  the  FXEM  Initializing  routines,  the  user  may  set  the  disposition 
of  FXEM  error  61  to  “continue"  rather  than  “abort.”  No  explicit  restart  features 
have  been  Incorporated  into  the  WFSM.  Such  procedures  would  not  be  feasible  In  an 
iterative,  space- stepped  solution. 


I 
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Appendix  A 


WFSM  FLOW  CHART 


O^'TOPT N?X 
\  Part  1. 
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Divide  accumulators  by 
total  distance  SACC 

VARBR  =  VACC  /  SACC 
VGBR  =  VGACC  /  SACC 
VCBR  =  VCACC  /  SACC 


Compute  distance-weighted 
route-mean  wind  factor 


=  VGBR  - 


/VCBR2+YARBR 

2*"ASPEED 


\ 

) 


WGARBR 


Appendix  B 


WFSM  PROGRAM  LISTING 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 
1  1 
12 

1  3 

14 

15 

16 

17 

18 
19 
2C 

2  1 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 
5  1 

52 

53 

54 

55 

56 

57 

58 

59 

60 


CENRWND  ENROUTE  WIND/R.  C.  WHITON/27  FEB  1979 
C 

SUBROUTINF  ENRWNO  (FRMLAT.  FRMING,  TOLAT.  TOLNG,  JULDAT, 
1  FCHRS.  IAIT.  IOPTN,  ASPEED,  GsPEED ) 


GMT. 


C* 

C* 

c* 

c* 

c* 

c* 

c* 


PROGRAM  ID- 
MET  ANALYST- 
SYS  aNALYST- 
PROGRAMMER- 

CREATED  ON- 


ENRWND 

MAJ  ROGER  C. 

MaJ  ROGER  C. 


WHITON,  USAFETAC/DNS .  EXT  5412 
WHITON,  USaFETaC/DNS,  EXT  5412 


MaJ  ROGER  C.  WHITON,  USaFETaC/DNS,  EXT  5412 


27  FEB  1979 


PROJECT- 


192301 


C* 

C*  DESCRIPTION- 
C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 


THIS  SUBROUTINE  SUBPROGRAM  CALCULATES  CLIMATOLOGI¬ 
CAL  WlNO  FACTORS  BY  SAWYER'S  EQUIVALENT  HEADWIND 
TECHNIQUE  FOR  REASONABLY  ARBITRARY  GREAT  CIRCLE 
ROUTES  AT  SPECIFIED  CONSTANT  ALTITUDES  FOR  ANY  OF 
FOUR  SEASONS  OF  THE  YEAR  IN  ANY  OF  FOUR  MODES  OR 
WIND  OPTIONS  (CALM  WIND,  SIMULATED  WIND,  MEAN 
WIND,  AND  90X  WORST  WIND).  CLIMB  AND  DESCENT 
WINDS  ARE  NOT  PROVIDED  FOR.  TEMPERATURES  AND 
AIRCRAFT  PERFORMANCE/FUEL  CONSUMPTION  ARE  NOT  PRO¬ 
VIDED  FOR. 


C* 


C* 

c* 

c» 

c« 

c» 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c« 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 


THE  SUBPROGRAM  NAVIGATES  A  SIMULATED  AIPCRAFT  OVER 
A  GREAT  CIRCLE  ROUTE  FROM  POINT  'A'  (FRMLAT,  FRMLNG) 
TO  POINT  'B'  (TOLAT.  TOLNG)  IN  A  GLOBAL  (LATITUDE/ 
LONGITUDE  ORIENTED)  GRID  SYSTEM  WHOSE  RESOLUTION  IS 
15  DEGREES  LATITUDE  OY  30  DEGREES  LONGITUDE.  THE 
SUBROUTINE  DETERMINES  WHAT  GRID  SECTOR  THE  AIRCRAFT 
IS  IN  AND  THE  LENGTH  OF  THE  AIRCRAFT'S  PATH  THROUGH 
THE  SECTOR.  THEN  THE  ROUTINE  CONSULTS  THF  WIND  DATA 
BASE  FOR  THE  PARTICULAR  GRID  SECTOR.  FLIGHT  LEVEL 
AND  SEASON  TO  OBTAIN  SECTOR-MEAN  WIND  DIRECTION, 
SPEED  AND  VARIANCE.  CROSS-TRACK  AND  ALONG-TRACK 
COMPONENTS  OF  THE  WIND  ARE  COMPUTED.  THESE  ARE 
WEIGHTED  BY  THE  LENGTH  OF  THE  AIRCRAFT'S  TRACK 
THROUGH  THE  SECTOR.  AT  THE  END  OF  THE  SIMULATED 
FLIGHT,  A  DISTANCE-WEIGHTED,  ROUTF-MEAN  CROSS-TRACK 
COMPONENT,  ALONG-TRACK  COMPONENT  AND  VARIANCE  IS 
COMPUTEC.  FROM  THESE,  A  ROUTE-MEAN  WIND  FACTOR  IS 
COMPUTED  AND  ADDED  TO  THE  AIR  SPEFD  TO  OBTAIN 
GROUND  SPEED.  L A T I TUDE /LONG  I TUOE  OF  POINTS  'A' 

AND  '8,'  THE  ALTITUDE  OF  THE  AIRCRAFT.  THF  JULIAN 
DATE  AND  GMT  OF  THE  FLIGHT.  AND  THE  AIRCRAFT'S 
AIR  SPEED  ARE  GIVEN.  THE  GROUND  SPEED  IS  CALCU¬ 
LATED  BASED  ON  WIND  FACTOR.  WIND  FACTOR  CALCULA¬ 
TIONS  CAN  BE  MADE  IN  ANY  OF  FOUR  MODES--CALM  WIND. 
SIMULATED  WIND.  90X  WORST  WIND.  AND  MEAN  WIND. 

IN  THE  CALM  WIND  CASE,  THE  ROUTE-MEAN  WIND  FACTOR 
IS  SIMPLY  SET  TO  ZERO,  WITH  NO  NAVIGATION  OF  THE 
AIRCRAFT  OR  ACCUMULATION  OF  D I S T ANCE -WE  I GHTE D 
WIND  COMPONENTS  BEING  NECESSARY.  THE  SIMULATED 
WIND  CASE  IS  RESERVED  FOR  FUTURE  DEVElOPMFnT,  IN 
WHICH  A  VARIABLE  WIND  MAY  BE  IMPLEMENTED  SUCH 
THAT  THE  WIND  FACTOR  NEED  NOT  Be  IDENTICALLY  THE 
SAME  EVERY  TIME  THE  ROUTE  IS  FLOWN.  PRESENTLY, 

THE  S I  MU).  A 1  ED  WIND  DEFAULTS  TO  THE  MEAN  WIND. 

IN  THE  MEAN  WIND  CASE,  THE  MEAN  WIND  IS  DFVELOPED 
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c* 
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80 
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c* 
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c* 
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c* 
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c* 
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c* 
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c* 

i  1  ? 

c* 
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c» 

1  1  4 

c* 

1  15 

c* 

l  16 

c* 
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c* 

1  10 

c» 

1  1  9 

c* 
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c* 

1  2  1 

c* 

1  22 

c* 

1  2  3 

c» 

124 

c* 

125 

c* 

1  26 

c* 

127 

c« 

l  2  8 

c* 

method- 


DIRECTLY  FROM  THE  MEAN  WIND  DATA  BASE.  IN  THE  90X 

WORST  CASE,  this  mean  wind  factor  is  statistically 

ADJUSTED  TO  A  VALUE  SUCH  THAT  ONLY  10X  OF  FLIGHTS 
OVER  THAT  ROUTE  WILL  EXPERIENCE  A  WORSE  WIND  FACTOR 
(  10X  RISK ) . 

THE  USER  DICTATES  WHETHER  HE  WANTS  CALM,  SIMULATED. 
MEAN  OR  90X  WORST  WINOS  BY  SPECIFYING  A  WIND  OPTION 
(IOPTN)  FROM  0  To  3  (SEE  TABLE  IN  PART  I  OF  SUBPRO¬ 
GRAM).  THE  USER  SPECIFIES  WHICH  OF  THE  FOUR  SEA¬ 
SONS  OF  THE  YEAR  HE  WANTS  WINDS  FOR  BY  STIPULATING 
A  JULIAN  DATE  (JULDAT)  FROM  1  TO  366.  THE  TABLE  FOR 
CONVERSION  OF  JULIAN  DATE  TO  SEASON  OF  THF  YEAR  IS 
PROVIDED  IN  PART  III  OF  THE  SUBPROGRAM.  THE  USER 
SPECIFIES  THE  ALTITUDE  OF  THF  AIRCRAFT  BY  SETTING 
THE  FLAG  I  ALT  TO  1  FOR  25.000  FT  AND  TO  2  FOR 
35.000  FT.  THE  USER  SPECIFIES  HIS  AIRSPEED  IN 
KNOTS  IN  THE  INPUT  VARIABLE  ASPEED  AND  RECEIVES 
CALCULATED  GROUND  SPEED  IN  KNOTS  IN  THE  OUTPUT 
VARIABLE  GSPEED. 

ARCHITECTURE  OF  THE  SUBPROGRAM  IS  SUCH  THAT  PART  II 
HANDLES  THE  CALM  WIND  CASE.  PART  III  THE  MEAN  WIND. 

PART  IV  THE  90/C  WORST  WIND.  AND  PART  V  THE  SIMULATED 

WINO.  PART  I  ACCOMPLISHES  BRANCHING  TO  PARTS  II-V. 
PART  VI  HANDLES  WRAPUP  AND  NORMAL  TERMINATION.  PtRT 

vii  handles  abnormal  termination. 

THE  METHODS  USED  BREAK  NEATLY  INTO  TWO  CATEGORIES. 
(1)  GEODESY  AND  NAVIGATION,  AND  (2>  METEOROLOGY. 

Cl)  GEODESY  AND  NAVIGATION... 

THE  AIRCRAFT  IS  NAVIGATED  ALONG  A  GREAT  CIRCLE  ROUTE 
BETWEEN  POINTS  'A*  AND  '6'  BY  MOVING  FROM  A  POINT 
CALLED  THE  ' CURRE NT-PO I N T •  (XLATC,  XLNGC >  TO  A  POINT 
CALLED  THE  'NEXT-POINT*  (XLATnP.  XLNGNP).  MOST  OF 
THE  COMPUTATION  OF  THE  NEXT-POINT  IS  DONE  IN  GLOBAL 
(LATITUDE /LONGITUDE)  COORDINATES. 

FOR  ANY  GIVEN  CURRENT-POINT,  THE  NEXT-POINT  IS 
EITHER  THE  POINT  WHERE  THE  GREAT  CIRCLE  ROUTE  CROS¬ 
SES  A  LATITUOF  GRID  LINE  (0.  15,  30.  ...  DEG)  OR 

THE  POINT  WHERE  THE  ROUTE  CROSSES  A  LONGITUDE  GRID 
LINE  (0.  30.  60.  ...  DEG),  WHICHEVER  IS  CLOSER  TO 

THE  CURRENT-POINT. 

AT  THE  CURRENT-POINT,  THE  AIRCRAFT'S  HEADING  CAN  BE 
CALCULATED  VIA  THE  HEADING  FORMULA  (FUNCTION  HDG > . 
KNOWLEDGE  OF  THE  LOCATION  OF  THE  C URR E N T -P 0  I N T  AND 
THE  AIRCRAFT'S  HEADING  AT  THAT  POINT  DICTATES  THE 
NEXT  LATITUDE  CROSSING  AND  THE  NEXT  LONGITUDE  CROS¬ 
SING.  CANDIDATE  #1  FOR  NEXT-POINT  IS  THE  POINT  OF 
NF  XT  LATITUOE  CROSSING  XLATNX  ( ANO  ASSOCIATED  LON¬ 
GITUDE  XLNdAS).  CANDIDATE  *2  IS  THE  POINT  OF  next 
LONGITUDE  CROSSING  XL  NGN  X  (AND  ASSOCIATED  LATITUDE 
XLATAS) . 

CALCULATION  OF  LONGITUDE  XLNGaS  ASSOCIATED  WITH  A 
GIVEN  LATITUDE  XLATNX  IS  DONE  ITERATIVELY  BY  SUB¬ 
ROUTINE  GRTCIR,  .HICH  SOLVES  THE  EQUATION  OF  A 
GREAT  CIRCLE.  NEWTON'S  ITERATIVE  METHOD  WITH  AN 
EXACT  SPHERICAL  TRIANGLE  FIRST  GUESS  IS  USED. 
virtually  all  of  the  complexity  OF  SUBROUTINES 

GRTCIR  ANO  FNRWN'-'  IS  ATTRIBUTABLE  IN  ONF  WAY  OR 
ANOTHER  TO  THIS  PROBLEM  OF  SOLVING  THE  EQUATION  OF 
A  GREAT  CIRCLE  FOR  LONGITUDF.  GIVEN  THE  LATITUDE. 
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calculation  of  latitude  xlatas  associated  with  a 

GIVEN  LONGITUDE  XLNGNX  IS  DONE  DETERMINISTICALLY 
8 Y  SUBROUTINE  GRTCIft.  NO  DIFFICULTY  IS  ENCOUNTERED 
HERE. 

CANDIDATES  FOR  NEXT-POINT  ARE  THEN  AS  FOLLOWS... 

CANDIDATE  a  1  ......  (XLATNX,  XLNGAS) 

CANDIDATE  V 2  .  (XLATAS.  XLNGNX) 

WITH  CANDIDATES  1  AND  2  IDENTIFIED,  SUBROUTINE 
DISTAN  IS  THEN  USED  TO  CALCULATE  GREAT  CIRCLE  DIS¬ 
TANCES  AS  FOLLOWS... 

51  ......  GCD  FROM  CURRENT-POINT  TO  CAND  «1 

5 2  .  GCD  FROM  CURRENT-POINT  TO  C AND  W2 

THE  CANDIDATE  WITH  THE  SHORTER  S-DISTANCE  IS  SE¬ 
LECTED  AS  THE  nEXT-PQInT. 

THE  CURRENT-POINT  CAN  THEN  OE  THOUGHT  OF  AS  THE 
POINT  WHERE  THE  AIRCRAFT  ENTERS  THE  GRID  SECTOR. 

AND  THE  NEXT-POINT  AS  THE  POINT  WHERE  THE  AIRCRAFT 
DEPARTS  THE  SECTOR.  THE  DISTANCE  SNP  BETWEEN 
CURRENT-POINT  AND  NEXT-POINT  IS  THE  DISTANCE  BY 
WHICH  THE  WIND  COMPONENTS  AND  VARIANCE  ARE  WEIGHTED 
IN  COMPUTING  THE  ROUTE-MEAN  WIND  FACTOR. 

2.  METEOROLOGY... 

WITH  CURRENT-POINT  AND  NEXT-POINT  ESTABLISHED.  AN 
APPROXIMATE  MIDPOINT  (XLATM,  XLnGM)  BETWEEN  THEM 
IS  ESTABLISHED.  THE  NUMBER  NSEC  INDICATING  WHICH 
OF  THE  108  GRID  SECTORS  THROUGH  THE  AIRCRAFT  IS 
PASSING  THROUGH  IS  THEN  CALCULATED.  BASED  ON 
(XLATM,  XLNGM).  SEE  DESCRIPTION  OF  THE  GRID  SYS¬ 
TEM  BELOW.  SINCE  FOUR  SEASONS  AND  TWO  FLIGHT 
LEVELS,  AS  WEIL  AS  108  GRID  SECTORS,  ARE  PROVIDED 
FOR.  A  DATUM-NUMBER  NRND  IS  OBTAINED  FROM  THE  GRID 
SECTOR  NUMBER  NS  C  (1*106).  THE  ALTITUDE  INDEX  IalT 
(1-2)  AND  THE  SEASON  INDEX  ISEASN  (1-4)  AS 
FOLLOWS.  .  . 

NRND  =  (ISEASN-l)  *  (TOTAL-NUMBER  ALTITUDES)  * 
(TOTAL-NUMBER-GRID-SFCTORS) 

+  (  I  ALT - 1  >  *  (TOTAL-NUMBER-GRID-SECTORs) 

♦  NSEC 

WITH  THE  DATUM-NUMBER  NRND  CALCULATED.  THF  DATA  BASE 
(SEE  DESCRIPTION  BELOW)  IS  CONSULTED  TO  OBTAIN  SEC¬ 
TOR-MEAN  WIND  DIRECTION  DIRM  (RADIANS,  IN  BETA-ANGLE 

form,  as  described  below),  speed  spdm  (knots)  and 

VARIANCE  VARM  ( KN0TS**2  )  . . . 

OIRM  =  DIR(NRND) 

SPDM  =  SPD(NRND) 

VARM  =  VAR(NRND) 

THE  AIRCRAFT  HEADING  AT  (XLATM,  XLNGM)  IS  COMPUTED 

by  means  of  Function  hog.  the  heading  hdgm  is  com¬ 
puted  in  beta-angle  format  (SFE  BFLOW)  and  in  that 

FORMAT  IS  REPRESENTED  BY  THE  ANGLE  ALPHA. 

ALPHA,  THE  WIND  DIRECTION  DIRM,  AND  THE  WIND  SPEED 
SPDM  ARE  USED  TO  CALCULATE  THE  ALONG-GROUND-TRACK 
WIND  COMPONENT  VGM  AND  THF  CROSS-GROUND-TRACK  WIND 
COMPONENT  VCM... 
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c* 

213 

c* 

2  14 

c* 

215 

c* 

2  16 

c* 

217 

c» 

216 

c* 

2  19 

c* 

220 

c» 

22  1 

c* 

222 

c* 

223 

c* 

224 

c* 

225 

c* 

226 

c* 

227 

c* 

220 

c* 

229 

c* 

230 

c* 

231 

c* 

232 

c* 

233 

c* 

234 

c* 

235 

c* 

236 

c* 

237 

c* 

230 

c* 

239 

c* 

240 

c* 

24  1 

c* 

242 

c* 

243 

c* 

244 

c* 

<45 

c* 

246 

c* 

24  7 

c* 

240 

c* 

<49 

c* 

<50 

c* 

251 

c* 

252 

c* 

<53 

c* 

<54 

C  4 

<55 

c* 

256 

c* 

<57 

C  ♦ 

<5R 

C* 

VGM  *  SPOM  *  COS(A8S<DIRM  - 
VCM  =  SPDM  *  SIN(  D I RM  - 

A  time-averaged  «ind  factor  for  the  present  leg  is 
then  computed  using  sawyer • S  EQUIVALENT  headwind 

FORMULA. . . 

WEAR  =  VGM  -  <(VCM**2  ♦  VARM/2.0)  / 

(2.0  *  ASPEED)) 

BUT  NO  FURTHER  USE  IS  MADE  OF  WOAR. 

THE  QUANTITIES  VGM.  VCM  AND  VARM  ARE  MULTIPLIED  BY 
THE  SECTOR-  OR  LFG-LENGTH  SNP.  AND  THF  PRODUCT  IS 
ACCUMULATED  IN  VGACC,  VCACC.  AND  VACC. 

THEN  THE  NEXT-POINT  BECOMES  THE  CURRENT-POINT,  A  NEW 
NEXT-POINT  IS  CALCULATED.  AND  THE  PROCESS  IS  RE¬ 
PEATED  FOR  THE  NEXT  LEG  OF  THE  SIMULATED  FLIGHT. 

AFTER  THE  LAST  LEG  OF  THE  FLIGHT  HAS  BEEN  PROCESSED. 
THE  ACCUMULATORS  VGACC.  VCACC  AND  V ACC  ARE  DIVIDED 
BY  THE  TOTAL  ROUTE  LENGTH.  THE  RESULTING  DISTANCE- 
WEIGHTED  ALONG-TRACK,  CROSS-TRACK  AND  VARIANCE  QUAN¬ 
TITIES  ARE  THEN  JSED  IN  THE  SAWYER  EQUIVALENT  HEAD¬ 
WIND  FORMULA  GIVEN  ABOVE  TO  OBTAIN  THE  DISTANCE- 
WEIGHTED  ROUTE-MEAN  WIND  FACTOR  WBARBR. 

IF  THE  90  X  WORST  WIND  FACTOR  HAS  BEEN  SELECTED* 
WBARBR  IS  STATISTICALLY  ADJUSTED  TO  THE  90X  WORST 
VALUE  BY  SAWYER'S  TECHNIQUE.  IMPLEMENTED  IN  PART  IV 
OF  SUBROUTINE  ENRWND. 

AS  PART  OF  THE  WTAPUP  ACTIONS  IN  PART  VI  OF  SUBROU¬ 
TINE  EnRWnO.  the  ROUTE-MEAN  wind  FACTOR  WBARBR  IS 
ADDED  TO  THE  AIR  SPEED  ASPEED  TO  GET  A  GROUND  SPEED 
GSPEEO  FOR  RETURN  TO  THE  PROGRAM  CALLING  SUBROUTINE 
ENRWNO. 

THE  METEOROLOGICAL  TECHNIQUE  USED  IS  REFERRED  TO  AS 
SAWYER'S  EQUIVALENT  HEADWIND  TECHNIQUE  AND  IS  DOCU¬ 
MENTED  IN  AWS-TR-77-267 ,  GUIDE  FOR  APPLIED  CLIMA¬ 
TOLOGY. 

BETA-ANGLE  FORMAT... 

THE  WIND  DIRECTION  DD  CAN  BE  EXPRESSED  IN  THE 
FORM  OF  A  BETA-ANGLE  B... 

I  * 

I  * 

I  * 

I  * 

I  * 

I  * 

i  *  e 


THE  BETA  ANGLE  DESCRIBES  THE  DIRECTION  TOWARD 
WHICH  THE  WIND  1 5.  BLOWING  IN  TERMS  OF  RADIANS 
COUNTERCLOCKWISE  FROM  EASTWARD.  CONVERSIONS 
ARE  AS  FOLLOWS... 


?4 


•W"  1  . .  I  w  II  mm  I  - 


259 

C* 

260 

C* 

DIRECTION  BETA-ANGLE  BETA-ANGLE 

261 

C* 

TOWARD  WHICH  (DEGREES)  (RADIANS) 

262 

C* 

263 

c* 

300.0  -30.0  -0.524 

264 

c* 

270.0  0.0  0.000 

265 

c* 

240.0  30.0  0.524 

266 

c* 

180.0  90.0  1.571 

267 

c* 

268 

c* 

GRID 

269 

c* 

system- 

THE  GRID  SYSTEM  IS  IN  GLOBAL  ( L AT  I TUDE /LOnG I TUDE > 

270 

c* 

COORDINATES.  WITH  A  RESOLUTION  OF  (15  DEG  LAT  X 

271 

c* 

30  DEG  LON).  THE  GRID  SYSTEM  EXTENDS  FROM  75  DEG 

272 

c* 

N  TO  60  DEG  S  (9  ROWS)  AND  GLOBALLY  IN  LONGITUDE 

273 

c* 

(12  COLUMNS).  HENCE,  THERE  ARE  108  GRID  SECTORS. 

274 

c* 

COLUMN  1  EXTENDS  FROM  30  DEG  W  TO  0  DEG  (GREEN¬ 

275 

c* 

WICH),  COLUMN  2  FROM  0  DEG  TO  30  DEG  E,  ETC.  ROW 

276 

c* 

1  EXTENDS  FROM  75  DEG  N  TO  60  DEG  N.  ROW  2  FROM 

277 

c* 

60  DEG  N  TO  45  DtG  N.  ETC.  SUCH  A  GRID  SYSTEM  IS 

276 

c* 

TOO  COARSE  FOR  CALCULATION  OF  OPERATIONAL  WIND 

279 

c* 

FACTORS  BUT  IS  SUITABLE  FOR  PURPOSES  SUCH  AS 

260 

c* 

SIMULATION. 

28  1 

c* 

282 

c* 

data 

283 

c* 

base- 

THE  SIOP  WINDS  WERE  USED  TO  CREATE  THE  DATA  BASE 

284 

c* 

EMPLOYED  BY  SUBROUTINE  ENRWND.  SIOP  WINDS  HAVE 

285 

c* 

A  (5  DEG  X  5  DEG)  RESOLUTION.  ALL  SIOP  WINDS 

286 

c* 

WITHIN  EACH  OF  THE  PRESENT  (15  DEG  X  30  DEG)  GRID 

287 

c* 

SYSTEM  WERE  EQUALLY  WEIGHTED  AnD  AVERAGED  TO 

286 

c» 

OBTAIN  A  SINGLE  WIND  for  EACH  OF  THE  108  GRID 

2  b  9 

c* 

SECTORS.  TWO  FLIGHT  LEVELS  (25,000  FT  AND 

290 

c* 

35,000  FT)  AND  FOUR  SEASONS.  THERE  ARE  THUS  A 

291 

c* 

TOTAL  OF  864  WINDS  IN  THE  DATA  BASE.  EACH  HAS  A 

292 

o 

WIND  DIRECTION  DIR(NRND)  (RADIANS.  BETA-ANGLE 

293 

c» 

FORMAT).  SPEED  SPD(NRND)  (KNOTS),  AND  VARIANCE 

294 

c* 

VAR(NRND)  ( KN0TS**2  > .  THESE  WINDS  ARE  STORED 

296 

c* 

SEQUENTIALLY  IN  THE  FOLLOWING  ORDER... 

296 

c* 

297 

c* 

WINTER  SEASON 

298 

c  ♦ 

FL  250 

299 

c* 

108  SECTORS 

300 

c* 

FL  350 

301 

c* 

108  SECTORS 

302 

c* 

SPRING  SEASON 

303 

c* 

FL  250 

304 

c* 

108  SECTORS 

305 

c* 

FL  350 

306 

c* 

106  SECTORS 

307 

c* 

•  •  •  • 

306 

c* 

ETC. 

309 

c* 

•  •  •  • 

310 

c* 

31  1 

c* 

THE  WINDS  ARF  In  THE  ARRAYS  DIRC870),  SPD(870) 

312 

c* 

AND  VARIB70).  THESE  ARRAYS  ARE  MEMBERS  OF  THE 

313 

c* 

LABELLED  COMMON  BLOCK  WEA... 

314 

c* 

315 

c* 

COMMON  /WEA/  DIRC870).  SPP(870>.  VAR(870> 

3  16 

c* 

317 

c* 

AN  INITIALIZING  ROUTINE  IS  NEEDED  TO  READ  THE 

3  1  8 

c« 

WIND  DATA  FROM  AN  864-RECORO  FILE  ON  DISK  OR 

3  1  9 

c* 

TAPE  TO  THE  WIND  ARRAYS.  THIS  INITIALIZATION 

320 

c* 

MUST  BE  EXECUTED  BEFORE  FIRST  CALL  TO  SUBROUTINE 

32  1 

c* 

ENRWND. 

322 

c* 

323 

c* 

limitations- 

PRESENTLY,  THE  WIND  FACTOR  SIMULATION  MODFL  CAN¬ 

324 

c* 

NOT  PRODUCE  A  'SIMULATED*  WIND  HAVING  A  DAY-TO- 

325 

c* 

DAY  VARIABILITY.  HENCE.  WITHIN  ANY  GIVEN  SEASON 

326 

c* 

ANO  FOR  ANY  GIVEN  ALTITUDE.  A  PARTICULAR  ROUTE 

25 


327 

C* 

WILL  ALWAYS  EXPERIENCE  THE  SANE  WIND  FACTOR,  RE¬ 

328 

C* 

GARDLESS  OF  THE  PASSAGE  OF  TIME.  REOUESTS  FOR 

329 

C* 

THE  SIMULATED  WIND  WILL  DEFAULT  TO  THE  MEAN 

330 

C* 

WIND. 

331 

C* 

332 

C* 

FURTHERMORE,  THE  PRESENT  MODEL  IS  INCAPABLE  OF 

333 

C* 

PRODUCING  A  'FORECAST •  WIND  FACTOR.  REQUESTS 

334 

C* 

FOR  A  FORECAST  WIND  WILL  DEFAULT  TO  THE  MEAN 

335 

C* 

WIND. 

336 

C* 

337 

C* 

TEMPERATURE.  AIRf-RAFT  PERFORMANCE,  FUEL  CONSUMP¬ 

336 

c* 

TION.  AND  OTHER  FACTORS  CONSIDERED  BY  TYPICAL 

339 

c* 

FLIGHT  PLANNING  «ODELS  ARE  NOT  INCLUDED  IN  THE 

340 

c* 

PRESENT  MODEL,  WHICH  DEALS  ONLY  WITH  WIND. 

34  1 

c* 

342 

c* 

OPERATIONS  NORTH  OF  75  DEG  N  OR  SOUTH  OF  60  0E6 

343 

c* 

S  ARE  NOT  PERMITTED. 

344 

c* 

345 

c* 

THE  COARSENESS  OF  THE  PRESENT  15  DEG  LATITUDE  BY 

346 

c* 

30  DEG  LONGITUDE  GRID  SYSTEM  PRECLUDES  USE  OF 

347 

c* 

THE  MODEL  FOR  CALCULATION  OF  OPERATIONALLY  REALIS¬ 

348 

c* 

TIC  WIND  FACTORS  UNLESS  CHANGES  ARE  MADE.  RATHER 

349 

c  * 

SIMPLE  MODIFICATIONS  TO  THE  EXISTING  CATA  BASE 

350 

c* 

AND  GRID  SYSTEM,  ALONG  WITH  CORRESPONDING  CHANGES 

351 

c* 

TO  THE  SOFTWARE,  CAN  ELIMINATE  THIS  LIMITATION. 

352 

c* 

353 

c* 

SPECIFICALLY.  THE  FOLLOWING  RESTRICTIONS  MUST  BE 

354 

c* 

ADHERED  TO... 

355 

C  ♦ 

356 

c* 

A.  NO  ROUTE  SHOULD  BE  FLOWN  DIRECTLY  OVER 

357 

c* 

EITHER  POLE. 

356 

c* 

359 

c* 

B.  NEITHER  POLE  SHOULD  SERVE  AS  A  POINT  'A' 

300 

c* 

(TAKEOFF)  OR  A  POINT  •!)'  (LANDING). 

361 

c* 

362 

c* 

c.  circumferential  or  round  robin  flights  in 

363 

c* 

WHICH  POINT  'A*  AND  '8*  COINCIDE  WILL  BE 

364 

c* 

aborted,  break  such  flights  into  smaller 

365 

c* 

SEGMENTS. 

366 

c* 

367 

c« 

D.  SEMl-CIRC'JMFERENTIAL  flights  in  WHICH  POINT 

368 

c* 

•B'  IS  EXACTLY  OPPOSITE  POINT  'A'  THROUGH 

369 

c« 

THE  CENTER  OF  THE  EARTH  WILL  ALSO  ABORT. 

370 

c» 

break  such  flights  into  smaller  segments. 

371 

c* 

372 

c* 

e.  routes  or  segments  of  routes  flown  directly 

373 

c* 

north  or  directly  south  along  a  longitude 

374 

c* 

line  will  also  abort. 

375 

c« 

376 

c* 

F.  AT  PRESENT.  THE  MODEL  ALLOWS  FLIGHTS  ONLY 

377 

c* 

AT  ALTITUDES  25,000  FT  ANO  35,000  FT. 

378 

c« 

379 

c* 

g.  all  departure  and  destination  points  must 

380 

c« 

LIE  within  the  latitudinal  range  From  75 

381 

c* 

DEG  N  TO  60  DEG  S. 

382 

c* 

383 

c* 

RESTRICTIONS  F  AND  G  WERE  Imposed  IN  ORDER  TO  RE¬ 

384 

c* 

DUCE  THE  CORE  STORAGE  REQUIREMENTS  OF  THE  MODEL. 

3b5 

c* 

386 

c* 

REFERENCES- 

1.  WHITON,  R.  C.,  AND  P.  l.  HERoo.  1980:  WIND  FACTOR 

387 

c* 

SIMULATION  MODEL:  MODEL  DESCRIPTION,  USAFETAC/TN-80/001 

3  d  8 

c« 

389 

c* 

390 

c* 

2.  HEROD,  p.  l . ,  and  r.  C.  lhiton.  1980:  WIND  FACTOR 

391 

c* 

SIMULATION  MODEL:  USER'S  MANUAL,  USAFETAC/TN-80/002 . 

392 

c* 

26 
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393 

c* 

394 

c« 

INPUT- 

FRMLAT  a  LATITUDE  OF  POINT  *A*  (TAKEOFF)  IN 

395 

c* 

DECIMAL  DEGREES.  NEGATIVE  SOUTH 

396 

c* 

FRMLNG  a  LONGITUDE  OF  POINT  'A*  (TAKEOFF)  IN 

397 

c« 

DECIMAL  DEGREES,  NEGATIVE  EAST 

39b 

c* 

TOLAT  a  LATITUDE  OF  POINT  'B*  (LANDING)  IN 

399 

c* 

DECIMAL  OEGREES,  NEGATIVE  SOUTH 

400 

c* 

T OLNG  *  LONGITUDE  OF  POINT  *B*  (LANDING)  IN 

401 

c* 

OECIMAL  DEGREES.  NEGATIVE  EAST 

402 

c* 

JULDAT  *  JULIAN  BASE  DATE  OF  WIND  FACTOR  RE¬ 

403 

c* 

QUEST.  INTEGER.  VALUES  1  -  366 

404 

c* 

GMT  *  GREENWICH  BASE  TIME  OF  WIND  FACTOR 

405 

c* 

REQUEST.  VALUE  0.0  -  24.0  HRS 

406 

c* 

FCHRS  *  FORECAST  HOURS  AHEAD.  VALUE  GREATER 

407 

c* 

THAN  OR  EQUAL  TO  0.0  HRS 

40b 

c* 

IALT  a  ALTITUDE  INDEX.  1  FOR  25,000  FT  AND 

409 

c* 

2  FOR  35.000  FT 

4  10 

c* 

IOPTN  a  WIND  OPTION  OR  MODE.  INTEGER.  0  FOR 

411 

c* 

CALM  WIND,  2  FOR  905C  WORST  WIND.  AND 

412 

c* 

3  FOR  MEAN  WIND.  IN  THE  FUTURE.  1 

413 

c* 

WILL  BE  IMPLEMENTED  AS  THE  SIMULATED 

4  1  4 

c* 

WIND  BUT  PRESENTLY  DEFAULTS  TO  THE 

4  15 

c* 

MEAN  WIND 

4  1  6 

c* 

ASPEED  a  ROUTE-MEAN  AIRCRAFT  AIR  SPEED  IN  KNOTS 

4  17 

c* 

4  1b 

c* 

0UTPUT- 

GSPEEO  a  ROUTE-MEAN  AIRCRAFT  GROUND  SPEED  IN 

4  1  9 

c* 

KNOTS 

420 

c* 

421 

c* 

system  sub¬ 

422 

c* 

programs 

423 

c* 

USED- 

SORT.  SIN.  COS.  FXEM 

424 

c* 

425 

c* 

USER  SUB¬ 

426 

c* 

PROGRAMS 

427 

c* 

USED- 

GRTCIR.  DISTAN.  SPHGLO,  HDG.  ORLAT.  8RLNG 

4  2  6 

c* 

429 

c* 

ESTIMATED 

430 

c* 

CPU  TIME- 

A  LONG  route  CONSISTING  OF  8  LEGS  AND  REQUIRING 

431 

c* 

AN  EQUATOR  CROSSING  TAKES  0.1  -  0.5  SEC  CPU  TIME 

432 

c* 

ON  A  HONEYWELL  6080  COMPUTER. 

433 

c* 

434 

c* 

STORAGE 

435 

c* 

rEquirements- 

PROCEDURE  PLUS  DATA  OCCUPY  1130  WORDS  OF  CORE 

4  36 

c* 

STORAGE.  COMMON  BLOCK  /WEA/  ADDS  2610  WORDS. 

437 

c* 

TOTAL  CORE  is  37*0  WORDS  for  this  SUBROUTINE. 

438 

c« 

439 

c« 

PROGRAM 

4  4  U 

c* 

updates- 

none 

44  1 

c* 

44? 

C***************************************************************«***| 

443 

c 

444 

COMMON  /WEA/ 

0IR(870),  SPDI07O),  V  AR ( 870 ) 

445 

c 

COMMON  /DBG/ 

NDE  BUG 

446 

c 

447 

c 

DESCRIPTION  of  SPHERICAL  GRID*  latitude  GRID  LINES  ARE  SPACED 

4  4  P 

c 

every  oeltla 

DEGREES  OF  LATITUDE.  LONGITUDE  GRID  LINES  ARE 

449 

c 

SPACED  EVERY 

DELTLO  DEGREES  OF  LONGITUDE. 

450 

c 

451 

DATA  DELTLA/15.0/.  DELTL0/3O.0/ 

452 

c 

453 

c 

angle  conversion  factors,  multiplier  dtor  converts  degrees  to 

454 

c 

radians,  rtod  CONVERTS  radians  to  degrees. 

455 

c 

456 

DATA  OTOR/O. 017*5329/,  RTOD/57 . 295 780/ 

457 

c 

458 

c 

LEGMAX  IS  maximum  PERMITTED  NUMBER  of  legs  IN  ROUTE,  used  to 

459 

c 

detect  runaway  routes. 

460 

c 

27 
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************************************* 

*  PART  I  * 

*  BRANCH  * 

A************************************ 

BRANCH  TO  APPROPRIATE  PART  OF  SUBROUTINE.  BASED  ON  VALUE  OF  IOPTN 
FLAG. 


FUNCTION 

CALM  WIND 
SIMULATED  WIND 
90X  WORST  WIND 
MEAN  WIND 


SUBROUTINE  PARTS 

PART  II 
PART  V 

PARTS  III  «  IV 
PART  III 


NOTE  THAT  PREPARATION  OF  THE  90V  WORST  WIND  REQUIRES  HAVING  PRE¬ 
VIOUSLY  ESTABLISHED  THE  MEAN  WInD. 

NOTE  THAT  IN  THIS  VERSION  OF  THE  SUBROUTINE.  THE  SIMULATED  WIND 
OPTION  IS  NOT  OPERATIONAL.  CALLS  FOR  THE  SIMULATED  WIND  WILL  DE¬ 
FAULT  TO  THE  MEAN  WIND. 

IF  < I OPTN  .LT.  0  .OR.  IOPTN  . GT •  3)  GO  TO  7020 
IOPTNP  =  IOPTN  ♦  1 

GO  TO  (2000.  5000,  3000.  3000).  IOPTNP 


PART  II 
CALM  WIND 


IN  THE  CALM  WIND  CASE,  THE  WIND  FACTOR  IS  ZERO.  AND  GROUND  SPEED 
IS  EQUAL  TO  A I p  SPEED  In  KNOTS.  RQUTE-A VERAGED  WInD  FACTOR  WBARBR 
IS  ZERO. 

WBARBR  a  0.0 

GO  TO  WRAPUP  AND  TERMINATION. 

GO  TO  6000 


PART  III 
MEAN  WIND 


INITIALIZE  ACCUMULATORS  FOR  DISTANCE  (SACC).  VARIANCE  (VACC). 
ALONG-TRACK  WIND  COMPONENT  (VGACC),  AND  CROSS-TRACK  WIND 
COMPONENT  (VCACC). 

SACC  =  0.0 
VACC  =  0.0 
VGACC  «  o.O 
VCACC  =  0.0 

INITIALIZE  LAST-LEG  FLAG  LASTLG.  WHEN  THF  CURRENT  LEG  nF  THE 
FLIGHT  IS  THE  LAST  LEG,  LAsTLG  a  1. 

LASTLG  =  0 

TEST  POSITION  OF  POINT  'A'  ( F Rml AT ,  F RMLNG )  AND  POINT  'B*  (TOLAT, 

tolnG)  for  legal  panGE  of  values  in  oegrefs. 


*>2  6 

IF 

(FRMLAT 

.LT. 

-90 . 0 

•  OR* 

FRMLAT 

.  GT  . 

90.0  ) 

GO 

TT 

7030 

527 

IF 

( TOLAT 

•  LT. 

-90.0 

•  Or* 

TOLAT 

•  G  r  • 

90.0) 

GO 

TO 

7030 

52R 

IF 

( frmlng 

•  LT. 

-180. 0 

•  OR. 

FRMLNG 

•  gt  • 

180.0) 

GO 

TC 

7030 

f 

I 

I 


529 

IF  (TOLNG 

.LT. 

-180.0 

•OR.  TOLNG  .GT.  180. 

0)  GO  TO  7030 

530 

C 

531 

C 

initialize 

INTERMEDIATE 

point  coordinates. 

532 

C 

533 

xlatc  ■  0 

•  0 

534 

XLNGC  =  0 

.0 

535 

THETAC  «  0 

•  0 

536 

PHIC  =  0 

•  0 

537 

C 

538 

c 

interpret 

THE  JULIAN  DATE  JULDAT  IN  TERMS  OF 

SEASON  OF  THE  YEAR 

539 

c 

Indicated 

By  the  season  Index  iseasn.  give 

an  error  condition 

540 

c 

FCR  ILLEGAL  JULDAT  VALUES. 

54  1 

c 

542 

c 

JUL I  an 

DATES 

SEASON  SEASON  INDEX 

543 

c 

JULDAT 

OF  YcAR 

iseasn 

544 

c 

545 

c 

1  - 

59 

W  I  NTC  R 

1 

546 

c 

60  - 

151 

SPRING 

2 

547 

c 

152  - 

243 

SUMMER 

3 

546 

c 

244  - 

334 

FALL 

4 

549 

c 

335  - 

366 

WINTER 

1 

550 

c 

551 

IF  (JULDAT 

.GT. 

366)  GO 

TO  7000 

552 

IF  (JULDAT 

.  G  E  . 

335)  GO 

TO  306b 

553 

IF  (JULDAT 

.  GE. 

244)  GO 

TO  3  0  5  j 

554 

IF  (JULDAT 

.GE. 

152)  GO 

TO  305b 

555 

IF  (JULDAT 

.GE  . 

60  )  GO 

TO  3060 

556 

IF  (JULOAT 

.GE. 

1  )  GO 

TO  3065 

557 

GO  TO  7000 

556 

3050 

ISEASN  =  4 

559 

GO  TO  3070 

560 

3055 

ISEASN  =  3 

56  1 

GO  TO  3070 

562 

3060 

ISEASN  «  2 

563 

GO  TO  3070 

564 

3065 

ISEASN  =  1 

565 

GO  TO  3070 

566 

C 

567 

C 

ChECK  ALTITUDE 

OPTION  IALT  FOR  LEGAL  RANGE. 

5  6  8 

c 

569 

c 

I  ALT 

AIRCRAFT  ALTITUDE 

570 

c 

571 

c 

1 

25*000  FT 

572 

c 

2 

35*000  FT 

073 

c 

574  3070  IF  (IALT  .LT.  1  • OP •  IAlT  .GT.  2)  GO  TO  7010 

57b  C 

576  C  CALCULATE  TOTAL  GREAT  CIRCLE  DISTANCE  STOT  IN  NAUTICAL  MILES  BE- 

577  C  TWEEN  POINTS  *A«  AND  «B.* 

b  7  B  C 


! 


579 
5b0 
58  1 
5  62 
583 

58  A 
585 
bfa  6 

587 

588 

589 

590 

59  1 

592 

593 

594 

595 

596 


5  T  0  T  =  0.0 

CALL  DISTAN  (FRMLAT ,  FRMLNG ,  TOLAT ,  TOL  NG  *  STOT) 

C 

c  initialize  great  circle  subroutine,  note  that  grtcir  is 

C  CALLED  TWICE--ONCE  FOR  ICONV  =  1  AND  ONCE  FOR  ICONV  =  2. 

C 

ICONV  =  1 

IGUFSS  =  0 
XLNGG  a  0.0 
XEAST  =  0.0 
XNORTH  =  0.0 

xlatc  =  o.o 

XLNGAP  =  0.0 
SAPRIM  =  0.0 

J075  CALL  GRTCIR  (ICONV.  FRMLAT,  FRMLNG,  TOLAT.  TOLNG .  XLATC.  XLNGC  , 
&  THETAA.  PHIA.  THETAB.  PHIB.  QTX.  OTY.  0T2,  IGUESS,  XLNGG. 

S  XEAST,  XNORTH,  XLNGAP,  SAPRIM,  XLATC.  THETAC,  PHIC) 

IF  (ICONV  .EQ.  2)  GO  TO  3080 


29 


4 


1 


597 

598 

599 

600 
60  1 
602 

603 

604 

605 

606 

607 

608 

609 

610 
611 
612 

613 

614 
6  15 
616 

617 

618 

619 

620 
62  1 
622 

623 

624 

625 

626 

627 

628 

629 

630 

631 

632 

633 

634 

635 

636 

637 

638 
6  39 

640 

641 

642 

643 


I  CON V  *  2 
GO  TO  3075 

SET  THE  LOCATION  ( XL ATC . XLNGC  )  OF  THE  CURRENT-POINT  TO  POINT 
•A*  (FRMLAT.FRMLNGl.  THIS  ROUTINE  OPERATES  WITH  TWO  PRINCIPAL 
POINTERS.  NAMELY  THE  CURRENT -PO I  NT  ( XLATC . XLNGC >  AND  THE  NExT- 
PO I  NT  (XLATNP.XLNGNP).  COMPUTATIONS  IN  THIS  SUBROUTINE  ARE 
PRIMARILY  IN  TERMS  OF  LATITUDE  AND  LONGITUDE.  SOME  LIMITED 
USE  OF  COLATITUDE  (THETA)  AND  AZIMUTH  (PHI)  IS  MADE  BY 
SUBROUTINE  GRTCIR.  SO  THESE  QUANTITIES  (THETAC.  PHIC)  ARE 
MAINTAINED  FOR  THE  CURRENT-POINT  (XLATC,  XLNGC). 

3080  XLATC  =  FRMLAT 
XLNGC  *  FRMLNG 
ICONVX  =  1 

CALL  SPHGLO  (ICONVX,  XLATC,  XLnGC,  THETAC,  PHIC) 

INITIALIZE  THE  LEG  COUNTER  LEGNO.  WHICH  IS  UNITY  FOR  THE  FIRST 

leg  of  the  flight. 


c 

C8000 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


legno  =  i 

IF  (NDEBuG  ,GT.  0)  WRITE  (6,8030) 

FORMAT  (1H0.  T  8 ■  ' CURR  PT*.  T21,  *NExT  PT*.  T34.  •DISTANCE'.  T44 , 

•SEC-1,  T50.  •  ACF  T  •  ,  T55.  • - WIND - ‘/IX. 

•LG’,  T  7 ,  •  lAT •  , 

T14,  'LNG1,  T  2  0 ,  'LAT*.  T27.  • LNG ' ,  T33,  'NAUT  MILES'.  T44, 
'TOR'.  T51,  ’HOG*.  T58.  *OIR'.  T64,  • SPD ' .  T70,  'VAR',  T74. 
•FACTOR • // ) 

BEGIN  NAVIGATION  LOOP. 

BASED  ON  A  GREAT  CIRCLE  BETWEEN  THE  CURRENT-POINT  AND  POINT  'B.' 
CALCULATE  THE  AIRCRAFT'S  HEADING  HDGC  AT  THE  CURRENT-POINT. 
HEADING  IS  THE  DIRECTION  TOWARD  WHICH  IN  DEGREES,  0.0  FOR  NORTH¬ 
WARD.  90.0  FOR  EASTWARD.  180. u  FOR  SOUTHWARD.  AND  270.0  FOR 
WESTWARD. 


3100  HDGC  =  HDG  (XLATC,  XLNGC.  TOLAT.  TOLNG) 

BASED  ON  THE  HEADING  HDGC  AT  T HE  CURRENT-POINT,  DETFRMlNE  THE 
EASTWARD  DIRECTION  INDICATOR  XE  AST  (POSITIVE  FOR  EASTWARD 
COMPONENT  OF  AIRCRAFT  MOTION)  AND  NORTHWARD  DIRECTION  INDI¬ 
CATOR  XNORTH  (POSITIVE  FOR  NORTHWARD  COMPONENT  OF  AIRCRAFT 
MOTION).  SPECIAL  CASES  OF  CARDINAL  HEADINGS  ARE  HANDLED  BY 
STATEMENTS  3110  THROUGH  3137.  WHILE  GENERAL  CASES  ARE  HANDLED 
BY  3140  THROUGH  3167. 

.EQ.  360.0)  GO  TO  3120 


64* 

3110 

IF  (HDGC  .60. 

0.0 

•  OR# 

HDGC 

645 

IF  (HDGC  .Eg. 

90.0  ) 

GO 

TO 

3125 

646 

IF  (HDGC  .EQ. 

180.0) 

GO 

TO 

3130 

647 

IF  (HDGC  .FO. 

270.0 ) 

GO 

TO 

3135 

646 

GO  TO  3140 

649 

3120 

XE AST  =  0.0 

650 

XNORTH  =  l.o 

651 

GO  TO  3170 

652 

3  125 

XEAST  =  1.0 

653 

XNORTH  =  0.0 

654 

GO  TO  3170 

655 

3130 

XEAST  =  0.0 

656 

XNORTH  =-l.o 

657 

GO  TO  3170 

656 

3135 

XEAST  *-1.0 

659 

XNORTH  =  0.0 

660 

3137 

GO  TO  3170 

661 

3140 

IF  (HDGC  .LT. 

90.0) 

GO 

TO 

3150 

662 

IF  (HDGC  .LT. 

180.0) 

GO 

TO 

3155 

663 

IF  (HDGC  .LT. 

270.0) 

GO 

TO 

3160 

664 

GO  TO  3165 

30 


665 

666 

667 

668 

669 

670 

671 

672 

673 

674 

675 

676 

677 
673 

679 

680 
66  1 
682 

683 

684 
665 
686 
687 
668 

689 

690 

691 

692 

693 

694 

695 
o96 

697 

698 

699 

700 

701 

702 

703 

704 

705 

706 

707 
706 

709 

710 

711 

712 

713 

714 

715 

716 

717 
7  1  8 

719 

720 
72  1 

722 

723 

724 

725 

726 

727 

728 

729 

730 

731 

732 


3150  XEAST  =  1.0 

XNORTH  *  1.0 

GO  TO  3170 
3155  XE AST  *  1.0 

XNORTH  =  -1.0 
GO  TO  3170 
3160  XE AST  =  -1.0 
XNORTH  =  -1.0 
GO  TO  3170 
3165  XE AST  =  -1.0 
3167  XNORTH  =  1.0 

3170  CONTINUE 

C  IF  (NDEBUG  ,GT.  11  WRITE  (6,8010)  XEAST ■  XNORTH 

C8010  FORMAT  (1H  ,  T 4 ,  • XE A S T . XNORTH  =  •,  T36.  2F10.4) 

C 

C  BASED  ON  EASTWARD  DIRECTION  INDICATOR,  CALCULATE  THE  POINT  OF  NEXT 

C  LONGITUDE  CROSSING  (XLNGNX)  AND  ASSOCIATED  LATITUDE  ( XL AT  AS )  DE- 

C  TERM  I  NED  FROM  THE  EQUATION  OF  A  GREAT  CIRCLE.  THIS  IS  CANDIDATE 

C  POINT  »1.  STATEMENTS  3200-3299.  DIFFERENT  PROCEDURE  IS  USED. 

C  DEPENDING  ON  WHETHER  THIS  IS  THE  FIRST  LEG  OF  THE  FLIGHT  OR 

C  WHETHER  (ON  SUBSEQUENT  LEGS)  THE  CURRENT-POINT  WAS  CANDIDATE  »1 

C  ( NC  AND  =  1)  OR  CANDIDATE  #2  ( NC  AND  =  2). 

C 

3200  IF  (LEGNO  .GT.  1)  GO  TO  3240 
C 

C  PROCEDURE  FOR  FIRST  LEG  OF  FLliHT.  DETERMINE  BRACKETING  LONGI- 

C  TUOE  VALUES.  MOVE  EAST  (DECREASING  LONGITUDE)  IF  XEAST  .GT. 

C  0.0.  SET  XLNGNX. 

C 

3210  CALL  8RLNG  (XLN6C,  XLNGLO,  XLNGHI) 

C  IF  (NDEBUG  .GT.  1)  WRITE  (6,8020)  XLNGLO.  XLNGHI 

C8020  FORMAT  <1H  .  TA.  • XLNGLO , XLNGH I  =  •,  T36.  2F10.A) 

IF  (XEAST)  3220,3220,3230 
3220  XLNGNX  =  XLNGHI 
Go  TO  3290 

3230  XLNGNX  *  XLNGLO 
GO  TO  3290 
C 

C  PROCEDURE  FOR  SUCCEEDING  LEGS  DEPENDS  ON  WHETHER  CANDIDATE  «1  OR 

C  CANDIDATE  #2  WAS  CHOSEN  LAST  TIME. 

C 

3240  GO  TO  (3250,3270),  NCAND 
C 

C  IF  CURRENT-POINT  WAS  A  CANDIDATE  81  (EXACT  LONGITUDE),  THEN  NEXT 

C  LONG  I TUOE  CROSSING  IS  THAT  LONGITUDE  PLUS  OR  MINUS  OELTA-LONGI- 

C  TUDE  (PLUS  FOR  WESTWARD,  MINUS  FOR  EASTWARD). 

C 

3250  XLNGNX  =  -(XEAST  *  DEL  TL  0  >  +  XlNGC 
C 

C  ADJUST  LONGITUDE  FOR  DATELINE  FOLD. 

C 

IF  (XLNGNX  .GT.  180.0)  GO  TD  J260 
IF  (XLNGNX  .LT.  -180.0)  GO  TO  3265 
GO  TO  3290 

3260  XLNGNX  =  XLNGNX  -  360.0 
GO  TO  3290 

3265  XLNGNX  =  XLNGNX  ♦  360.0 
GO  TO  3290 
C 

C  IF  CURRENT-POINT  WAS  A  CANDIDATE  82  (INEXACT  LONGITUDE).  THEN 

C  NEXT  LONGITUDE  CROSSING  MUST  8fc  FOUND  BY  THE  BRACKETING  METHOD 

C  USED  FOR  LEGNO  =  1. 

C 

3270  GO  TO  3210 
C 

c  set  associated  latitude  (xlatas). 

c 

3290  ICONV  «  3 


31 


733 

CAUL  GRTCIR  (ICONV.  FRML  AT  •  FRMLNG.  TOLAT.  T0LN6.  XLATAS, 

XLNGNX. 

734 

i 

THETAA.  PHIA»  THE  T  A8 •  PH1B,  GTX,  OTY.  0T2.  IGUFSS.  XLNGG. 

735 

& 

XEAST.  XNORTM,  XLNGAP,  SAPRIM.  XLATC.  ThETAC.  PHIX) 

736 

c 

IF  ( NOEBUG  .GT.  1)  WRITE  <6, 80301  XLATAS.  XLNGNX 

737 

CS030 

FORMAT  (1H  ,  T4.  •  XLATAS  .XLnGNX  =  T36.  2F10.4) 

738 

c 

739 

c 

BASED  ON  NORTHWARD  OIRECTIOn  INDICATOR.  CALCULATE  THE  POINT  OF 

740 

c 

NEXT  LATITUOE  CROSSING  (XLATNX)  AND  ASSOCIATED  LONGITUDE 

741 

c 

(XLNGAS)  DETERMINED  FROM  THE  EQUATION  OF  A  GREAT  CIRCLE. 

THIS 

742 

c 

IS  CANDIDATE  POINT  »2.  STATEMENTS  3300-3399.  DIFFERENT 

PRO- 

743 

c 

CEOURE  IS  USED.  DEPENDING  WHETHER  THIS  IS  FIRST  LEG  OF  FLIGHT 

744 

c 

OR  WHETHER  (ON  SUBSEQUENT  LEGS)  THE  CURRE NT-PO I  NT  WAS  A 

745 

c 

CANDIDATE  hi  (NCAND  =  1>  OR  A  CANDIDATE  #2  ( NC  AND  =  2). 

746 

c 

747 

c 

FIRST  SET  COUNTER  OF  • GR  TC I R •  CONVERGENCE  FAILURES  TO  ZERO 

• 

748 

c 

THIS  WILL  LATER  BE  USED  TO  REVERSE  THE  LaTITUDF  SEARCH. 

749 

c 

750 

3300 

NF  AIL  =  0 

751 

C 

75? 

C 

DISCRIMINATE  FIRST  FROM  SUBSEQUENT  LEGS. 

753 

C 

754 

3310 

IF  (LEGNO  .GT.  1)  GO  TO  3340 

755 

C 

756 

C 

PROCEDURE  FOR  FIRST  LEG.  DETERMINE  BRACKETING  LATITUDE  VALUES. 

757 

C 

move  north  (Increasing  latitude)  if  xnorth  .gt.  o.o.  set 

758 

C 

xlatnx. 

759 

C 

760 

3315 

CALL  BRLAT  (XLATC.  XLATLO.  XLATHI) 

761 

c 

IF  (N0E8UG  .GT.  1)  WRITE  (6,6040)  XLATLO.  XLATHI 

76? 

C  8  0  4  0 

FORMAT  (1H  .  T  4 ,  'XLATLO. XLATHI  =  T36.  2F10.4) 

763 

IF  (XNORTH)  3320.3320.3330 

764 

3320 

XLATNX  =  XLATLO 

765 

GO  TO  3390 

766 

3330 

XLATNX  =  XLATHI 

767 

GO  TO  3390 

7  o  8 

C 

769 

C 

PROCEDURE  For  succeeding  legs  depends  ON  WHETHER  CANDIDATE 

m  or 

770 

C 

CANDIDATE  *2  WAS  CHOSEN  LAST  TIME. 

771 

C 

77? 

3340 

GO  TO  (3350.3370),  NCANO 

773 

C 

774 

C 

IF  CURRENT-POINT  WAS  A  CANDIDATE  «l  (INEXACT  LATITUDE),  THEN  NEXT 

775 

c 

LATITUOE  CROSSING  MUST  BE  FOUND  BT  THE  BRACKETING  METHOD 

USED  FOR 

776 

c 

LEGNO  =  1. 

777 

c 

778 

3350 

GO  TO  3315 

779 

C 

780 

c 

if  current-point  was  a  candidate  ta  (Exact  latitude),  then 

next 

781 

c 

latitude  crossing  is  that  latitude  plus  OR  minus  dflt a-latitude 

78? 

c 

(PLUS  FOR  northward,  MINUS  FOR  EASTWARD). 

78? 

c 

7P4 

3370 

xlatnx  =  (XNORTH  *  DELTLA)  ♦  XLATC 

785 

IF  (  A8S( XLATNX )  .LE.  0.001  )  XLATNX  =  0.0 

786 

C 

i 

767 

C 

ADJUST  LATITUDE  FOR  POLAR  SINGULARITY. 

78P 

c 

789 

IF  (XLATNX  .GT.  90.0)  GO  TO  3375 

j 

790 

IF  (XLATNX  .LT.  -90.0)  GO  TO  3380 

791 

GO  TO  3390 

1 

792 

3375 

XLATNX  =  180.0  -  xlatnx 

1 

793 

GO  TO  3390 

794 

3380 

XLATNX  =  -180.0  -  XLATNX 

1 

795 

GO  TO  3390 

796 

C 

797 

C 

ATTEMPT  TO  SET  ASSOCIATED  LONGITUDE  (XLNGAS). 

798 

C 

799 

3390 

ICONV  =  5 

bOO 

C 

IF  (NOEBUG  .GT.  1)  WRITE  (6,8050)  XLATNX 

c 

C8060 

C 

CS070 

C 

C8080 

c 

c 

c 

c 

c 

c 

c 


c 

C  8  090 


FORMAT  (1H  ,  TA.  • XL ATNX  =  T36.  F10.4) 

CALL  GRTCIR  (ICONV.  FRMLAT.  FRMLNG,  TOLAT.  tolng.  xlatnx, 
xlngas.  thetaa,  phia.  thetab.  phib.  gtx.  qty.  qtz. 
iguess.  xlngg,  xeast.  xnortm,  xlngap.  saprim.  xlatc,  ThETAX. 

PHIC  > 

RESCUE  PHIC  FOR  LATER  USE. 

iconvx  =  i 

CALL  SPHGLO  (ICONVX,  XLATC.  XLNGC,  THETAC.  PHIC) 

IF  (NOEBUG  .GT.  1)  WRITE  (6,8060)  ICONV 
FORMAT  (1H  .  T  4 ,  ‘ICONV  =  T36.  110) 

IF  (NDE8UG  .GT.  1)  WRITE  (6,8070)  XLNGG 

FORMAT  (1H  .  T  4 ,  ‘FIRST  GUESS  LONGITUDE  =  T36.  F10.4) 

IF  (NOEBUG  .GT.  1)  WRITE  (6,8080)  XL ATNX ,  XLNGAS 
FORMAT  (IH  ,  T  4 ,  • XL ATNX , XLNGAS  *  •,  T36.  2F10.4) 

IF  'GRTCIR'  SUBROUTINE  FAILS  TO  CONVERGE  IN  LONGITUDE,  INCREMENT 
FAILURE  COUNTER  BY  ONE,  REVERSE  LATITUDE  SEARCH  DIRECTION.  AND 
TRY  AGAIN  TO  DETERMINE  THE  NEXT  LATITUDE  CROSSING  (XLATNX)  AND 
ASSOCIATED  LONGITUDE.  OTHERWISE.  CONTINUE.  TWO  CONSECUTIVE 
FAILURES  PRODUCE  AN  ERROR  ABORT. 

IF  (ICONV  .NE.  11)  GO  TO  3400 
NFAIL  =  NFAIL  ♦  1 
IF  (NFAIL  .GT.  1)  GO  TO  7040 
X  NOR  T  h  =  -XNORTH 

IF  (NOEBUG  .GT.  1)  WRITE  (6,8090)  XNORTH 

FORMAT  (1H  ,  T  4 ,  'REVERSED  XNORTH  =  ••  T36.  F10.4) 

GO  TO  3310 

CALCULATE  GREAT  CIRCLE  DISTANCE  (GCD)  si  from  current-point 
(XLATC. XLNGC)  TO  CANDIDATE  Ml  ( XL  AT  AS . XLNGNX )  .  CALCULATE 
GCD  S2  FROM  CURRENT-POINT  TO  CANDIDATE  M2  (XLATNX, XLNGAS). 

CALL  D1STAN  (XLATC.  XLNGC.  XL  A  T  A  S ,  XLNGNX,  SI) 

CALL  0 1  ST  AN  (XLATC.  XLNGC.  XLATNX,  XLNGAS.  $2 ) 

WHICHEVER  OF  THE  TWO  POINTS,  CANDIDATES  1  OR  2,  IS  CLOSEST  TO  THE 
CURRENT-POINT  IS  SELECTED  AS  THE  NEXT-POINT  ( X L A T NP , X L NGNP ) 


1 

C 

along 

the  Path.  DECISION  IS 

made  on  the  basis 

8  4  2 

C 

DIDATE 

POINT  #1  AND  S  2 

FOR  CANDIDATE  POINT 

*2. 

843 

C 

844 

IF  <  S2 

- 

Si )  3420*3420* 

34  30 

84b 

c 

846 

c 

CASE  OF 

Si  • GE •  S  2  •  •  • 

CHOOSE 

CANDIDATE  * 2  • 

847 

84  b 

3420 

NCAND 

= 

2 

34  9 

X  L  A  T  N 

= 

XL  ATNP 

6  80 

XLNGN 

= 

XLNGNP 

o5  1 

XLATNP 

= 

XLATNX 

ab? 

XLNGNP 

= 

XLNGAS 

«53 

SNP 

= 

S? 

65  4 

GO  TO  3450 

8  5  5 

C 

656 

C 

case  of 

S  2  •  GE  .  si  ... 

CHOOSE 

CANDIDATE  HI. 

ob  7 

c 

6b8 

3430 

NCAnD 

X 

1 

859 

xlatn 

= 

XLATNP 

860 

XLNGN 

= 

XLNGNP 

86  1 

XLATNP 

S 

XL  ATAS 

86  2 

XLNGNP 

= 

XLNGNX 

fr  63 

SNP 

= 

S  1 

8  6  4 

C 

665 

c 

NOTE  •  •  • 

( XLATN* XLNGN  ) 

IS  CROSSING  POINT  BEFORE 

866 

c 

H67 

c 

CALCULATE  GREAT  circle 

D1ST ANC 

E  (GCD)  SNEXT 

IN  1 

FROM  POINT  'A*  (FRMLAT.FRMLNG)  TO  THE  NEXT-POINT  ( XLATNP , 


33 


»69 

Q70 

871 

872 
073 

874 

875 

876 

877 
078 

879 

880 
88  1 
682 
083 

884 

885 
686 
087 
008 

889 

890 

891 

892 
o9  3 
094 
095 
096 

897 

898 
39  9 

900 

901 
9C2 

903 

904 

905 

906 

907 

908 

909 
9  10 
9  1  1 
9  1  c 
9  1  3 
9  1  4 
9  1  5 
916 
9  1  7 
9  1  fc 
9  1  9 
920 
92  1 
92  2 

923 

924 
92  8 
92  6 

927 

928 

929 

930 

931 

932 

933 

934 

935 

936 


C  XLNGNP)  Just  selected.  IF  that  distance  <snext>  is  greater 

C  THAN  THE  total  GCD  (STOT)  BETWEEN  'A'  AND  •B.*  THEN  THE  NEXT- 

C  POINT  IS  BEYOND  THE  TERMINATION-POINT  'B.*  IN  SUCH  A  CASE. 

C  THE  COORDINATES  OF  »B*  MUST  BE  USED  AS  THE  LAST-POINT 

C  (XLATNP, XLNGNP).  DISTANCE  SNP  MUST  BE  ADJUSTED.  AND  LAST-LEG 

C  FLAG  MUST  BE  TURNED  ON. 

C 

3450  SnExT  =  0.0 

C  IF  (NDEBUG  ,GT.  1)  WRITE  (6.8100)  NC  AND .  XLATNP.  XLNGNP 

CfllOO  FORMAT  (1H  ,  T4,  • NCAND , XL ATNP , XLNGNP  =  *.  T36,  110,  2F10.4) 

CALL  DISTAN  (FRMLAT,  FRMLNG,  XLATNP,  XLNGNP,  SNEXT) 

IF  (STOT  -  SNEXT)  3460.3500.3500 
C 

C  PROCEDURE  FOR  SNEXT  .GT.  STOT. 

C 

3460  LASTLG  =  1 

XLATNP  =  TOLAT 
XLNGNP  =  TOLNG 

CALL  DISTAN  ( XL ATN •  XLNGN ,  TOLAT.  TOLNG.  SNP) 

C 

C  FIND  APPROXIMATE  MIDPOINT  BETWEEN  CURRE N T - PO I N T  ( XL ATC . XLNGC ) 

C  AND  NEXT-POINT  ( XL A TNP . XLNGNP ) .  THIS  MIDPOINT  IS  LATER  USED 

C  TO  DETERMINE  GRID  SECTOR  NUMBER  NSEC  AND  AS  POSITION  FOR 

C  COMPUTATION  OF  AIRCRAFT  HEADING  HDGM. 

C 

3500  XLATM  =  0.5  *  (XLATC  ♦  XLATNP) 

XLNGM  =  0.5  *  (XLNGC  ♦  XLNGNP) 

C 

C  FIND  ROW  NUMBER  NROW  (LATITUDE  COUNTER)  AND  COLUMN  NUmBFR  NCOL 

C  (LONGITUDE  COUNTER)  OF  MIDPOINT.  THEN  USE  NROW  AND  NCOL  TO 

C  COMPUTE  SECTOR  NUMBER  NSEC  IN  >'HICH  THE  MIDPOINT  (XLATm, 

C  XLNGM)  LIES.  IF  POINT  LIES  OUTSIDE  PLUS  OR  MINUS  75.0  DEGREES 

C  LATITUDE  LIMITS  OF  THE  GRID  SYSTEM,  BOGUS  THE  POINT  INTO  THE 

C  SYSTEM. 

C 

XlATDM  =  XLATM 

IF  (XLATDM  .GT.  75.0)  XLATDM  =  75.0 

IF  (XLATDM  .lT.  -60.0)  XLATDM  =  -60.0 
ANGLE  =  75.0  -  XL ATOM 
NROW  =  (ANGLE  /  DELTLA)  ♦  I 
ANGLE  =  -XLNGM  +  30.0 

IF  (XLNGM  .GT.  0.0)  ANGLE  =  390.0  -  XLNGM 
NCOL  =  (ANGLE  /  DELTLO )  ♦  I 
IF  (NCOL  .GT.  1?)  NCOL  =  NCOL  -  12 

NSEC  =  ((NROW  -  I)  *  12)  ♦  1  *  (NCOL  -  1) 

C 

C  CALCULATE  DATUM  NUMBER  NRND  BASED  ON  SECTOR  NUMBER  NSEC 

C  (1-106).  AIRCRAFT  ALTITUDE  IAlT  (1-2)  AND  SEASON  INDEX 

C  ISEASN  (1-4). 

c 

NRND  =  (ISEASN  -  1)  *  2  *  108  ♦  (IALT  -1)4  108  +  NSEC 

C 

c  calculate  aircraft  heading  hdg<  (direction  toward  which  in 

C  DEGREES)  AT  THE  MIDPOINT  *M.' 

c 

HDGM  =  HOG  (XLATM.  XLNGM .  TOLAT,  TOLNG) 

C 

C  CONVERT  HEADING  IN  OEGRFES  TO  ANGLE  ALPHA  IN  RADIANS.  ALPHA 

C  IS  THE  DIRECTION  TOWARD  WHICH  THE  AIRCRAFT'S  GROUND  TRACK  IS 

C  INCLINED.  IN  RADIANS  COUNTERCLOCKWISE  FROM  EASTWARD. 

C 

ALPHA  =  90.0  -  HOGM 

IF  (ABS(HDGM)  -  180.0)  35  fc  0 . 35  f>0 . 35  5  0 

3550  ALPHA  =  360.0  +  ALPHA 
3560  ALPHA  =  0  T  OR  *  ALPHA 
C 

C  RETRIEVE  OIRECTION  OIRm  AND  SPEED  SPDM  TOWARD  WHICH  THE  WIND  IS 

C  3L0WING  AT  THE  MIDPOINT  POSITION  (XLATM, XLNGM).  DIRECTION  IS 
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937 

938 

939 

940 
94  1 

942 

943 

944 

945 

946 
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960 
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962 

963 
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967 
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969 
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978 
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900 
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982 
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986 

987 
96b 

989 

990 
99  1 

992 

993 

994 

995 

996 

997 

998 

999 
1  000 
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1003 
1  004 


C  IN  RADIANS  COUNTERCLOCKWISE  F-<OM  EASTWARD  (E.G..  WIND  FROM  225 

C  DEGREES  HAS  DIRM  =  +PI/4  =  0.78540).  SPEED  IS  IN  KNOTS.  ALSO 

C  RETRIEVE  VECTOR  STANDARD  VARIANCE  VARM  (SQUARE  OF  VECTOR  STAN- 

C  DARD  DEVIATION)  IN  UNITS  OF  KN0TS**2. 

C 

DIRM  x  DIR(NRNO) 

SPDM  =  SPD(NRND) 

VARM  x  VAR(NRND) 

C 

C  COMPUTE  ALONG-GROUND-TRACK  WIND  COMPONENT  V6M  AND  CROSS-GROUND- 

C  TRACK  WIND  COMPONENT  VCM  IN  KNOTS. 

C 

GAMMA  x  DIRM  -  ALPHA 

GAMMAA  x  AS  S ( GAMMA ) 

VGM  x  SPDM  *  CDS ( GAMMAA ) 

VCM  x  SPOM  *  SIN1GAMMA  ) 

C 

C  COMPUTE  TIME-AVERAGED  WIND  FACTOR  W8AR  IN  KNOTS  FOR  THE 

C  PRESENT  LEG. 

C 

WBAR  =  VGM  -  (IVCMWVCM  +  (VARM/2.0) )  /  (2.0WASPEED)  ) 

C  IF  (NDEBU6  .GT.  0)  WRITE  (6,8110)  LEGNO.  XL  AT  C ,  XLNGC ,  XLATNP. 

C  XLNGNP .  SNP,  NSEC,  HDGM.  DIRM.  SPDM,  VARM,  WBAR 

CdllO  FORMAT  (1H  , 12 , 2( F6. 1 ,F7. 1 ) , T3 3 , F8 . 2 , T4 2 , 14, lx, 2F7. 3.F6. 1 , 

C  F6.0,  F7.1) 

C 

c  ACCUMULATE  total  distance  to  THE  ACCUMULATOR  SACC.  distance- 

C  WEIGHTED  total  VARIANCE  TO  V  A  C  C  ,  DISTANCE-WEIGHTED  TOTAL 

c  alonG-track  component  to  vgacc,  and  distance-weighted 

c  total  cross-track  component  to  vcacc. 

c 

SACC  =  SACC  +  SNP 

VACC  =  VACC  ♦  (SNP  *  VARM) 

VGACC  x  VGACC  ♦  (SNP  *  VGM  ) 

VCACC  =  VCACC  ♦  (SN>'  *  VCM  ) 

c 

c  advance  to  next-point,  old  next-point  i xlatnp. xlngnp >  becomes 

C  CURRENT-POINT  ( X L AT C . XL NGC ) .  UPDATE  (ThETAC,  PHIC)  COR- 

C  RESPONDING  TO  (XLATC,  XLNGC).  INCREMENT  LEG  NUMBER  COUNTER. 

C  GIVE  ERROR  CONDITION  FOR  RUNAWAY  ROUTE  (EXCESSIVE  NUMBFR  OF 

C  LEGS). 

C 

XLATC  =  XLATNP 
xlngc  x  xlngnp 
iconvx  =  i 

CALL  SPHGLO  (ICONVX,  XLATC.  XL  ,GC ,  THETAC.  PHIC) 

LEGNO  =  LEGNO  +  1 

IF  (LEGNO  .GT.  LEGMAX)  GO  TO  7,50 
C 

C  IF  CURRENT-POINT  ( XL A T C . X L N G C )  IS  POINT  *B’  ( T 0 L AT , T OL NG ) ,  THIS 

C  IS  THE  LAST  LEG  OF  THE  FLIGHT,  AND  SUBROUTINE  BRANCHES  TO  WRAP- 

C  UP.  OTHERWISE,  ROUTINE  GOES  BACK  FOR  THE  NEXT  LEG. 

C 

xlaoif  =  assixlatc  -  tolat> 

XLNDIF  x  AB3(XLNGC  -  TOLNG) 

IF  ( XL  AD  I F  .LT.  0.0001  .AND.  XLNDIF  .LT.  0.0001  )  LASTLG  *  1 
IF  (LASTLG  .EC.  0)  GO  TO  3100 
C 

C  COMPUTE  FINAL  route-averaged  wind  factor  WBARBR.  DIVIDE  DISTANCE- 

C  WEIGHTED  ACCUMULATORS  BY  TOTAL  DISTANCE  TO  GET  ROUTE-MEAN  VARI- 

C  ANCE  VAR6P.  ROUTE-MEAN  ALONG-TRACK  COMPONENT  VGBR ,  AND  ROUTE- 

C  MEAN  CROSS-TRACK  COMPONENT  VCBR.  THEN  CALCULATE  THE  ROUTE-MEAN 

C  WIND  FACTOR  WBARBR. 

C 

VARBR  *  VACC  /  SACC 
VGBR  X  VGACC  /  SACC 
VC8R  x  VCACC  /  SACC 

WBARBR  x  VGBR  -  ((VCBR4VCBR  +  (VARBR  /  2.0))  /  (2.0  *  ASPEED)) 
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c 

C8120 

C 

C 

c 

c 

c 

c 

c 

c 


IF  (NDE8UG  .GT.  0)  WRITE  (6,8120)  STOT,  VGBR.  VCBR,  VARBR,  WBARBR 
FORMAT  <1Mo.  'GREAT  CIRCLE  DISTANCE*.  T30.  F15.5//IX. 
•ROUTE-AVERAGE'.  TI6.  • AL0N6-TRACK • .  T30.  FIS. 5/ 

T16.  'CROSS-TRACK',  T30.  FI5.S/T16,  'VARIANCE',  T30.  FlS.S/ 

T 1 6 ,  'WIND  FACTOR*.  T30.  F15.S) 

IF  IOPTN  =  2.  BRANCH  TO  PART  IV  OF  THE  SUBROUTINE  TO  ADJUST  THE 
WIND  FACTOR  TO  THE  90*  WORST  CASE.  OTHERWISE,  BRANCH  TO  PART 
VI  FOR  WRAPUP  ANO  NORMAL  TERMINATION. 

GO  TO  (6000.  6000,  *000.  6000).  IOPTNP 

************************************* 

*  PART  IV  * 

*  90*  WORST  WIND  * 

WWW********************************** 

90*  WORST  WIND  FACTOR  CALCULATION  EMPLOYS  THE  SAWYER  METHOD 
DESCRIBED  In  AIR  WEATHER  SERVICE  TECHNICAL  REPORT  77-267, 

GUIDE  TO  APPLIED  CLIMATOLOGY.  PAGE  6-8.  THE  NEW  VARIABLES 
USED  IN  THIS  CALCULATION  ARE  AS  FOLLOWS... 

XKFACT  *  FACTOR  TO  CONVERT  MEAN  STANDARD  VECTOR  DEVIATION 
OF  winds  over  a  route  <STVDVW)  to  STANDARD  devia¬ 
tion  of  THE  MEAN  WIND  FACTOR  (STDVWF) 

STVDVW  =  MEAN  STANDARD  VECTOR  DEVIATION  OF  WINDS  OVER  A 
ROUTE 

STOVWF  =  STANDARD  DEVIATION  OF  THE  MEAN  WIND  FACTOR 
WBARBR  =  90*  WORST  WIND  FACTOR  (I.E.,  10*  RISK) 

CALCULATE  K-FACTOR  OF  SAWYFR.  THE  RELATION  USED  BELOW  IS  A 

geometric  Curve  fit  to  sawyer's  data  reported  in  awstr  77- 

267.  THE  GEOMETRIC  CURvE  FIT  WAS  CONSTRUCTED  USING  hEwLETT- 
PACKARO  65  STAT  RAC  #2.  PAGE  22.  THE  FUNCTION  GIVES 
K  (XKFACT)  AS  A  FUNCTION  OF  ROUTE  LENGTH  (STOT). 

XKFACT  =  0.70585066  *  ( 0 . 999030 56** S TO T ) 

CALCULATE  MEAN  STANDARD  VECTOR  DEVIATION  OF  WIND  OVER  ROUTE. 
STVDVW  =  SORT ( VARBR ) 

calculate  standard  deviation  of  the  mean  wind  factor. 

STCVWF  =  STVDVW  *  XKFACT 
CALCULATE  907.  WORST  WIND  FACTOR. 

WBARBR  =  WBARBR  -  (1.28  *  STDV.F) 

WRAPUP  AND  NORMAL  TERMINATION. 

GO  TO  6000 


PART  V 

S I MUL ATED  WIND 


simulated  wind  presently  defaults  to  the  mean  wind  case. 


GO  TO  3000 


♦  PART  VI  * 

*  WRAPUP  t  NORMAL  TERMINATION  * 

•••*••*•*•*•»***«***•*•****•********* 

V I  NO  factor  is  ground  speed  minus  air  speed,  use  this  relation 

TO  COMPUTE  GROUND  SPEED  (GSPEED)  IN  KNOTS  FROM  GIVEN  AIR  SPEED 
(ASPEED)  IN  KNOTS  AND  ROUTE “AVERAGED  HIND  FACTOR  (WBARBR)  IN 
KNOTS. 

GSPEED  =  ASPEED  ♦  WBARBR 
NORMAL  TERMINATION. 


*  PAnT  VII  * 

*  abnormal  termination  * 

************************************* 


GSPEED  »  ASPEED 

WRITE  <6.8130)  JULDAT.  GMT.  IALT,  FRMLAT.  FRMLN6,  T0LAT. 
TOLNG.  IOPTN 

format  (  1  HO .  '»***» InD  factor  FRROR  D I AGNO ST l C S* * ** • / l X . 
•JULIAN  DATE  ■  T 1 5 ,  110.  T35,  'GMT  =  ••  TA 1 .  F10.2.  T53. 

•ALT  OPTION  =  •.  T66.  I10/1X. 

•FROM-LAT  *  ••  FIO. 2.  T  30 .  • FRQM-LNG  s  *,  TA 1 .  F10.2/1X. 

•TO-LAT  =  *.  FIO. 2.  T32.  •  TO-LNG  «=  •.  Flo. 2.  T52. 

•WIND  OPTION  *  110//) 

CALL  FXEM  <61.  'SUBROUTINE  ENR  WND • . . ILLEGAL  JULIAN  DATE*.  10) 
RETURN 

GSPEEO  =  ASPEED 

WRITE  <6.8130)  JULDAT,  GMT.  IALT.  FRMLAT.  FRMLNG,  TOLAT. 
TOLNG.  IOPTN 

CALL  FXEM  <61.  'SUBROUTINE  ENRWND. . . ILLEGAL  ALTITUDE*.  9) 
RETURN 

GSPEED  =  ASPEED 

WRITE  <6.8130)  JULDAT.  GMT.  IALT,  FRMLAT.  FRMLNG,  TOLAT, 
TOLNG.  IOPTN 

CALL  FXEM  (61.  'SUBROUTINE  E NR» ND . . . I LL EG AL  WIND  OPTION',  10) 
RETURN 

GSPEED  =  ASPEED 

WRITE  <6,8130)  JULDAT,  GMT,  IALT,  FRMLAT,  FRMLNG,  TOLAT. 
TOLNG,  IOPTN 

CALL  FXEM  <61.  'SUBROUTINE  ENR . ND . . . I LL EG AL  LAT/LON'.  9) 
RETURN 

GSPEEO  =  ASPEED 

WRITE  (6,8130)  JULDAT.  GMT,  IALT,  FRMLAT,  FRMLNG,  TOLAT, 
TOLNG,  IOPTN 

CALL  FXEM  <61.  'SUBROUTINE  E NR WND . . . GR T C I R  FAILED  TWICE*.  10) 
GSPEED  =  ASPEED 

WRITE  (6,8130)  JULDAT.  GMT,  IALT,  FRMLAT,  FRMLNG.  TOLAT. 
TOLNG.  IOPTN 

CALL  FXEM  <61.  'SUBROUTINE  ENR WND ...  RUN A WAV  ROUTE*.  9) 

RE  TURN 

end 


1  CGRTCIR  EG  GRT  CIRCLE/R.  C.  WH1T0N/03  FEB  1979 

2  C 

3  SUBROUTINE  GRTCIR  (ICONV,  FRMLaT,  FRMLNG.  TOLAT,  TOLNG.  XLAT. 

4  &  XLNG •  THETAAt  PHIA,  THETAB.  PHIB.  OTX,  QTY.  QTZ.  IGUESS.  XLNGG. 

5  &  XEAST,  XNORTH,  XLNGAP.  SaPRIM,  XL  ATC  «  THETA.  PHI) 

6  C 

7  C******** *********************************** **************************** 


8 

C* 

9 

C* 

PROGRAM 

10- 

GRTCIR 

10 

C* 

MET  ANALYST- 

MAJ  ROGER  C.  WHITON,  USAFETAC/DNS .  EXT  5412 

1  1 

c* 

SYS  ANALYST- 

MAJ  ROGER  C.  WHITON,  USAFETAC/DNS.  EXT  5412 

12 

c* 

PROGRAMMER- 

MAJ  ROGER  C.  WHITON,  USAFETAC/DNS.  EXT  5412 

13 

c* 

14 

c* 

CREATED 

ON- 

05  FEB 

1979  PROJECT-  192301 

15 

c* 

16 

c* 

DESCRIPTION- 

THIS  SUBROUTINE  SUBPROGRAM  SOLVES  THE  EQUATION  OF 

1? 

c* 

THE  GREAT  CIRCLE  DEFINED  BY  POINTS  *  A '  (  ( FRML AT . 

18 

c* 

FRMLNG)  OR  ( THETAA.PHI A) )  AND  <B*  (( TOLAT , TOLNG ) 

19 

c* 

OR  ( THETAB. PHIB ) )  AT  AN  INTERMEDIATE  POINT 

20 

c* 

( <  XL AT 

.XLNG)  OR  ( THETA, PHI )) .  DEPENDING  ON  THE 

21 

c* 

VALUE 

OF  THE  FLAG  ICONV,  THE  PROGRAM  PERFORMS  ONE 

22 

c* 

QF  THE 

Following  functions... 

23 

c* 

24 

c* 

ICONV 

function 

25 

c* 

26 

c* 

1 

GIVEN  INPUT  LATITUDE  AnD  LONGITUDE  OF  POINTS 

27 

c* 

•A*  (FfiMLAT. FRMLNG)  AND  'B*  ( TOLAT . TOLNG ) . 

28 

c* 

CALCULATE  SPHERICAL  COORDINATES  1  THE  T  AA , 

29 

c* 

PHIA)  AND  (THETAB, PHIB)  OF  THE  POINTS  By 

30 

c* 

INVOKING  subroutine  sphglo.  then  compute 

31 

c* 

THE  CROSS  PRODUCT  VECTOR  QT  =  OA  X  QB  OF 

32 

c* 

•A*  AND  •B.' 

33 

c* 

34 

c* 

2 

INITIALIZE  THE  EXACT  SPHERICAL  TRIANGLE 

35 

c* 

METHOD  OF  OBTAINING  FIRST  GUESS  LONGITUDE. 

36 

c* 

FIRST,  DETERMINE  THE  LONGITUDE  XLnGAp  OF 

37 

c* 

THE  POINT  A-PRIME  AT  WHICH  THE  GREAT  CIR¬ 

38 

c* 

CLE  ROUTE  CROSSES  THE  EOUATOR.  THEN  ES¬ 

39 

c* 

TABLISH  THE  INTERIOR  ANGLE  A-PRIME  OF  THE 

40 

c* 

SPHERICAL  TRIANGLE  ( A-PRIME .C-PRIME ,B ) . 

41 

c* 

THE  SINE  OF  THAT  ANGLE  IS  SAPRIM. 

42 

c* 

43 

c* 

3 

GIVEN  INPUT  INTERMEDIATE  LONGITUDE  XLNG, 

44 

c* 

CALCULATE  INTERMEDIATE  PHI,  INTERMEDIATE 

45 

c* 

THETA  AND  INTERMEDIATE  LATITUDE  XL AT . 

46 

c* 

47 

c* 

4 

GIVEN  INPUT  INTERMEDIATE  PHI,  CALCULATE  IN¬ 

48 

c* 

TERMEDIATE  THETA  AND  INTERMEDIATE  LATITUDE 

49 

c« 

XL  AT  . 

50 

c* 

51 

c* 

5 

GIVEN  INPUT  INTERMEDIATE  LATITUDE  XL  AT , 

52 

c* 

CALCULATE  INTERMEDIATE  THETA,  INTERMEDIATE 

53 

c* 

PHI  AND  INTERMEDIATE  LONGITUDE  XLNG. 

54 

c* 

55 

c* 

6 

GIVEN  INPUT  INTERMEDIATE  THETA,  CALCULATE 

56 

c* 

INTERMEDIATE  PHI  AND  INTERMEDIATE  LONGI¬ 

57 

c* 

TUDE  XLNG. 

58 

c* 

59 

c* 

THIS  SUBROUTINE  MUST  BE  INVOKED  ONCE  FOR  FUNCTION 

60 

c* 

ICONV 

*  1  AND  ONCE  FOR  ICONV  =  2  BEFORE  BEING 

61 

c* 

USED  FOR  FUNCTION  ICONV  =  3  THROUGH  6.  THUS 

62 

c* 

THERE 

MUST  BE  'INITIAL*  CALLS  TO  THIS  ROUTINE 

63 

c* 

USING 

ICONV  =  1  AND  2  BEFORE  CALLS  USING 

64 

c* 

ICONV 

*  3  THROUGH  6. 

65 

c* 

66 

c« 

METHOD- 

THE  METHOO  USED  EMERGES  FROM  THE  VECTOR  CALCULUS 

67 

c* 

IN  A  RECTANGULAR  (X.Y.Z)  COORDINATE  SYSTEM.  QA  IS 

68 

c* 

THE  POSITION  VECTOR  OF  POINT  'A.'  WHILE  08  IS  THAT 

38 


OF  POINT  'B.'  The  position  vector  of  THE  INTER¬ 
MEDIATE  POINT  IS  0.  THE  PLANE  OF  THE  GREAT  CIRCLE 
IS  DEFINED  BY  THE  CROSS  PRODUCT  QT  *  QA  X  OB. 

WHICH  IS  DIRECTED  NORMAL  TO  THE  PLANE.  THE  EQUA¬ 
TION  OF  THE  GREAT  CIRCLE  IS  DEFINED  BY  THE  RELATION. 

Q  .  (QA  X  QB)  =  0 

WHERE  THE  PERIOD  INDICATES  A  DOT  PRODUCT.  THIS  IS 
EQUIVALENT  TO  SAYING  THAT  THE  COMPONENT  OF  Q  IN  THE 
QT-0 1 REC  T ION  IS  ZERO.  EXPANSION  OF  THIS  PRODUCT  IN 
RECTANGULAR  COORDINATES  LEADS  TO  THE  RELATION, 

QTX*SIN(THETA)*COS(PHI  )  ♦  QT Y * S I N ( THE T A ) * S I N ( PH  1 ) 

♦  QTZ*COS( THETA )  =  0 

WHERE  QTX ,  QTY  AND  QTZ  ARE  THE  X,  Y  AND  Z-COMPONENTS 
OF  THE  CROSS  PRODUCT  QT.  DIVIDED  BY  R**2.  THE 
SQUARE  OF  THE  RADIUS  OF  THE  EARTH.  SOLVING  THIS 
EQUATION  FOR  THETA,  THE  COLATITUDF.  LEADS  TO  A 
DIRECTLY  SOLVABLE  EQUATION, 

THETA  =  ATAN(-QTZ/(QTX*COS(PHI )  ♦  QT y* S I N < PH  I > > 

ON  THE  OTHER  HAND,  SOLVING  FOR  PHI  LEADS  TO  A 
TRANSCENDENTAL  EQUATION  OF  THE  FORM, 

Q  T  X*COS ( PH  I  >  ♦  QTy*SIN(PHI)  ♦  Q TZ /T AN ( THETA )  =  0 

THIS  LATTER  EOUATION  IS  SOLVED  BY  NFWTON'S  ITERA¬ 
TIVE  TECHNIQUE  FOR  NON-LINEAR  EQUATIONS.  THE 
CRITICAL  FIRST  GUESS  FOR  PHI  IS  PROVIDED  BY  ONE 
OF  THREE  METHODS.  DEPENDING  ON  USER  DESIRES  AND 
WHETHER  THE  POINT  A-PRIM  HAS  ALREADY  BEEN  ESTAB¬ 
LISHED.  ANY  ONE  OF  THE  THREE  FIRST  GUESS  METHODS 
CAN  BE  USED.  THE  METHODS  ARE... 

(1)  USER-SPECIFIED  FIRST  GUESS 

(2)  APPROXIMATELY  PLANAR  TRIANGLE 

(3)  EXACT  SPHERICAL  TRIANGLE. 

NORMALLY  THE  THIRD  METHOD  IS  USED.  BUT  METHOD  (2) 
MUST  BE  USED  AT  LEAST  ONCE  TO  BREAK  GROUND  FOR 
(3).  METHOD  (2)  IS  ONLY  APPROXIMATE  AND  LOSES 
ACCURACY  THE  LONGER  THE  DIS¬ 
TANCE  BETWEEN  *  A  *  AND  'B.'  THIS  CAN  BE  A  MAJOR 
LIMITING  FACTOR,  BECAUSE  NEWTON'S  TECHNIQUE  REQUIRES 
A  REASONABLY  CLOSE  FIRST  GUESS  FOR  PHI. 

ONS-  1.  THE  APPROXIMATELY  PLANAR  TRIANGLE 

METHOD  USED  AS  FIRST  GUESS  FOR  NEWTON'S  ITERATIVE 
TECHNIQUE  APPEARS  TO  BREAK  DOWN  WHEN  ABSOLUTE  VAL- 
OF  THE  DELTA-LONGITUDE  OR  DELTA-LATITUDE  BETWEEN 
POINTS  '  A '  AND  '*}'  EXCEEDS  ABOUT  40-50  DEGREES,  PAR¬ 
TICULARLY  WHEN  THE  TWO  POINTS  ARE  IN  DIFFERENT  QUAD¬ 
RANTS  OF  THE  GLOBE.  UNDER  THESE  CIRCUMSTANCES.  NEW¬ 
TON'S  METHOD  CONVERGES  TO  A  FALSE  ROOT  OR  FAILS  TO 
CONVERGE  ENTIRELY,  AND  IT  BECOMES  IMPOSSIBLE  TO 
SOLVE  THE  EQUATION  OF  A  GREAT  CIRCLE  FOR  THE  AZIMUTH 
ANGLE  PHI.  POINTS  'A'  AND  'B'  ALONG  THE  GREAT 
CIRCLE  MUST  THUS  BE  KEPT  REASONABLY  CLOSE  TOGETHER. 
IF  POINT  'B'  IS  TOO  FAR  REMOVED  FROM  'A.'  THEN  USE 
THIS  SUBROUTINE  DETERMINISTICALLY  (ICONV  =  I  AND  2 
FOLLOWED  BY  ICONV  *  3  OR  4 >  TO  CALCULATE  THETA 
( XL  A  T )  FOR  AN  INTERMEDIATE  PHI  (XLNG)  NOT  MORE  THAN 
30  DEGREES  FROM  PHIA.  THEN  USE  THAT  POINT  AS  THE 
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61 

62 

63 


67 


70 

71 

72 

73 

74 

75 

76 

77 

78 

79 
l  80 
181 
18? 
183 
1  84 
185 
1  86 
1  87 
188 
1  89 
1  90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

00 

01 

02 

03 

04 


C* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c» 

c* 

c* 

c* 

c* 

c* 

c* 

C*  INPUT- 
C* 

C* 

c* 

c* 

c* 

c* 

c 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c» 

c* 

c* 

c* 

c* 

c» 

c* 

C*  INPUT  i 
C*  OUTPUT- 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 


NEW  POINT  8  IN  ALL  FURTHER  GREAT 
TIONS.  THIS  CAN  BE  DONE  BECAUSE 
POINT  NEED  NOT  BE  BOUNDED  BY  'A* 


CIRCLE  COMPUTA- 
THE  •INTERMEDIATE' 
AND  •B.* 


2.  IN  THE  CASE  OF  A  GREAT  CIRCLE  THAT  IS  ALSO 
A  LONGITUDE  CIRCLE.  THE  IMPLICIT  FUNCTION 
F I  THETA • PH I >  =  0.0  DOES  NOT  RETURN  A  UNIQUE 
SOLUTION  IN  THETA  FOR  A  SPECIFIED  PHI.  IF  THE 
GREAT  CIRCLE  IS  ALONG  A  LONGITUDE.  ANY  VALUE 
OF  THETA  IS  A  SOLUTION  TO  THE  FUNCTION.  THIS 
SUBROUTINE  CANNOT  BE  USED  TO  CALCULATE  THETA 
FROM  PHI  IF  FRMLNG  =  TOLNG.  AN  ERROR  MESSAGE 
IS  DISPLAYED  ON  SYSOUT  AND  ERROR  FLAG  IC0NV*12 

is  set  if  this  is  attempted. 


frmlat 
FRMLNG 
TOLAT 
TOLNG 
THET AA 

PHIA 


THETAB 

PHIB 


LATITUDE  OF  POINT  'A'  (ORIGIN)  IN  DEG 

LONGITUDE  OF  POINT  *A'  (ORIGIN)  IN  DEG 

LATITUDE  OF  POINT  'B'  (DESTINATION)  IN  DEG 

LONGITUDE  OF  POINT  '8*  (DESTINATION)  IN  DEG 

COLATITUDE  ANGLE  IN  RADIANS  OF  POINT  *A,' 
0.0  TO  PI  RADIANS  SOUTHWARD  FROM  NORTH. 
POSITIVE 

AZIMUTH  ANGLE  IN  RADIANS  COUNTERCLOCKWISE 
FROM  GREENWICH  LOOKING  DOWNWARD  AT  THE 
NORTH  POLE.  PLUS  OR  MINUS  0.0  TO  PI  RADIAN 
(NEGATIVE  MEANS  CLOCKWISE) 

COLATITUDE  IN  RADIANS  OF  POINT  *8' 

azimuth  in  radians  of  point  *b' 


IGUESS  =  flag  CONTROLLING  WHETHER  FIRST  GUESS  IS 
SPECIFIED  FOR  LONGITUDE  In  VARIABLE 
XLNGG.  IF  IGUESS  =  1.  XLNGG  IS  INTER¬ 
PRETED  AS  FIRST  GUESS  LONGITUDE.  IF 
IGUESS  =  0.  XLNGG  Is  IGNORED  ON  INPUT 
HUT  IS  FILLED  WITH  THE  PROGRAM-GENERATED 
FIRST  S.tSS  VALUE  ON  OUTPUT. 


XE AST  =  EASTWARD  DIRECTION  INDICATOR,  POSITIVE 
FOR  AIRCRAFT  MOTION  TOWARD  THE  EAST. 


XNORTH  =  NORTHWARD  DIRECTION  INDICATOR,  POSITIVE 
FOR  AIRCRAFT  MOTION  TOWARD  THE  NORTH. 


XLATC  =  LATITUDE  OF  CURRENT-POINT,  I.E..  THE 

POINT  FROM  WHICH  THE  PROGRAM  IS  MARCHING 
IN  OBTAINING  XLAT,  THE  LATITUDE  OF  THE 
SO  CALLED  NEXT-POINT. 


ICONV  *  FLAG  CONTROLLING  FUNCTION  TO  BE  PERFORMED 
BY  THIS  SUBROUTINE.  SEE  TABLE  IN  •METHOD* 
ABOvF.  WHEN  RETURNED  AS  OUTPUT.  ICONV  =  1 
THROUGH  6  FOR  NORMAL  TERMINATION  OR  A  NUM¬ 
BER  GREATER  THAN  OR  EOUAt  TO  10  FOR  ABNOR¬ 
MAL  TERMINATION.  ABNORMAL  TERMINATION 
CODES  FOR  ICONV  ARE... 


ERRONEOUS  ICONV  SPECIFIED 
NEWTON'S  ITERATIVE  TECHNIQUE  FOR 
SOLUTION  OF  NON-LINEAR  EOUATTONS 
DID  NOT  CONVERGE  IN  NM \X  IT!  R  A - 


1  0 
1  1 


SYSTEM  SUB¬ 
PROGRAMS 
USED- 

USER  SUB¬ 
PROGRAMS 
USED- 

ESTIMATED 
CPU  TIME- 

STORAGE 

REQUIREMENTS 

PROGRAM 

UPDATES- 


TIONS 

12-  ATTEMPT  WAS  MADE  TO  CALCULATE 
THETA  FROM  PHI  FOR  A  FLIGHT 
ALONG  A  LONGITUDE  CIRCLE 

*  LATITUDE  OF  INTERMEDIATE  POINT  IN  DEG 
«  LONGITUDE  OF  INTERMEDIATE  POINT  IN  DEG 

*  X-COMPONENT  OF  CROSS  PRODUCT  QT  DIVIDED  BY 
R**2.  COMPUTED  IN  FUNCTION  ICONV  =  1  AND 
REUSED  IN  SUBSEQUENT  CALLS 

a  Y-COMPONENT  OF  ABOVE 

=  Z-COMPOmENT  OF  ABOVE 

=  FIRST  GUESS  LONGITUDE  FOR  ICONV  =  4 
OR  S.  IF  VALUE  IS  NOT  SPFCIFIED  AS 
INPUT.  THE  PROGRAM-GENERATED  FIRST 

guess  will  de  provided  as  output 


XLNGAP  =  LONGITUDE  IN  DEGREES  OF  THE  POINT  A-PRIME 
AT  WHICH  THE  GREAT  CIRCLE  CROSSES  THE 
EQUATOR 

SAPRIM  =  SINE  OF  INTERIOR  ANGLE  A-PRIME  OF  THE 
EXACT  SPHERICAL  TRIANGLE 

THETA  =  COLATITUDE  OF  INTERMEDIATE  POINT  IN  RADIANS 

PHI  =  AZIMUTH  OF  INTERMEDIATE  POINT  IN  RADIANS 

NOTE-  ALL  SOUTH  LATITUDES  AND  EAST  LONGITUDES  ARE 
REPRESENTED  AS  NEGATIVE  NUMBERS. 


SIN,  ARSIN.  COS,  ARCOS.  TAN,  ATAN,  FXEM 


SPHGLO,  HOG 


TO  BE  SUPPLIED. 


-  PROCEDURE  PLUS  DATA  OCCUPY  1042  WORDS  OF  STORAGE. 


EPSILON  IS  THE  CRITERION  FOR  CONVERGENCE  OF  NEWTON'S  ITERATIVE 
TECHNIQUE  AND  HAS  UNITS  OF  RADIANS.  PI  IS  THE  GEOMETRIC  PI. 
NMAX  IS  THE  MAXIMUM  NUMBER  OF  ITERATIONS  PERMITTED  FOR  THE  NEW¬ 
TON  technique,  ntrifg  is  a  flag  indicating  whfther  the 
DIRECTION  (TRIDIR)  OF  THE  SPHERICAL  TRIANGLE  IN  FIRST 
GUESS  METHOD  »3  HAS  BEEN  REVERSED  DUE  TO  AN  EQUATOR 
CROSSING.  NORMALLY  ZERO,  NTRIFG  =  l  INDICATES  A  REVERSAL. 

COMMON  /DBG/  NDF  BUG 

DIMENSION  APRMD<2>.  XLNGGG  <  2 ) ,  PH  I GG ( 2  > 

DATA  EPS/O. 0001/,  P I / 3 . 1 4 1 5927/ ,  NTRIFG/O/.  NMAX/7/ 


273  C  BRANCH  TO  APPROPRIATE  SECTION  OF  THE  SUBROUTINE  BASED  ON  THE  FUNC- 

274  C  TION  REQUESTED  VIA  THE  ICONV  FLAG.  THE  ERROR  CONDITION  ICONV*IO 

275  C  IS  RETURNFP  FOR  ILLEGAL  ICONV  VALUES.  IF  FLIGHT  IS  ALONG 

276  C  A  LONGITUOE  CIRCLE  AND  ICONV  REQUESTS  CALCULATION  OF  THETA. 

277  C  AN  ERROR  CONDITION  IC0NV=12  IS  SET  AND  AN  ERROR  MESSAGE  IS 

276  C  OISPLAYED  On  SYSOUT.  ICONVS  IS  A  SAVED  VALUE  OF  THE  ORIGINAL 

279  C  ICONV. 

280  C 


281 

ICONVS  *  ICONV 

282 

IF 

( ICONV  . GE • 

1 

.  ANO. 

ICONV 

,LE.  5) 

GO 

TO 

50 

283 

ICONV  =  10 

284 

GO 

TO  1000 

285 

50  IF 

(ICONV  . NE . 

2 

•  AND. 

ICONV 

.NE.  3) 

GO 

TO 

60 

286 

IF 

( ABS ( FRMLNG 

- 

TOLNG 1 

•  GT  . 

0.0001 ) 

GO 

TO 

60 

287 

288 

289 

290 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 

301 

302 

303 

304 

305 

306 

307 

308 

309 

310 
31  1 

312 

313 

314 


ICONV 


12 


GO  TO  1020 

60  GO  TO  (100,  200.  300,  400,  500,  600).  ICONV 


************************************* 

*  PART  I  * 

*  FUNCTION  ICONV  =  1  * 

*  CROSS  PRODUCT  COMPUTATION  * 


CONVERT  (LAT.LNG)  OF  *A*  AND  *R*  TO  ( THE T A A , PH  I A )  AND  (THETA8, 
PHI8  )  . 

100  ICONVX  =  1 

CALL  SPHGLO  (ICONVX.  FRMLAT.  frmlng.  THETaA.  PHIA) 
call  SPHGLO  (ICONVX,  TOLAT,  TOLNG.  TheTAB.  PhIB) 

PRELIMINARIES  FOR  COMPONENTS  OF  POSITION  VECTORS  QA  AND  OB. 

SINTHA  =  SIN(THETAA) 

COSTHA  =  CqS(THETAA) 

SINPHA  =  SINIPHIA  1 
COSPHA  =  COS (PHI  A  ) 

sinthb  *  sin(tmetab) 

COSTHB  =  COS(THETAB) 

SlNPHB  =  SIN(PHIB  ) 

CQSPHB  =  COS(PHIB  1 


Jib 

316 

31  7 

318 

319 

320 

32  1 

322 

323 

324 


C 

C  RECTANGULAR  COMPONENTS  OF  POSITION  VECTORS  OA  AND  OP,  DIVIDED 

C  BY  THE  RADIUS  OF  THE  EAPTH  R. 

C 

OAX  =  SINTHA  *  COSPHA 
GAY  =  SINTHA  *  SINPHA 
0  A  2  =  COSTHA 

GBX  =  SInTHB  *  COSPHB 
OB  Y  r  SINTHB  *  SlNPHB 
QBZ  =  COSTHB 


325 

C 

326 

C 

COMPUTE 

COMPONENTS  of 

the 

CROSS  PRODUCT  VECTOR 

QT  = 

QA  X 

327 

C 

DIVIDED  BY  THE  RADIUS  OF 

THF  EARTH  SQUARED,  R**2. 

3  2  B 

c 

329 

QTX  = 

Q  A  Y  *  QBZ 

QAZ 

*  09  Y 

330 

OT  Y  - 

QAZ  ♦  3BX 

OAX 

♦  QBZ 

331 

QTZ  = 

Q  A  X  *  QB  Y 

QA  Y 

*  QB  X 

332 

c 

333 

c 

INITIALIZE  TRIANGLE  FLAG, 

USED  TO  REVERSE  THE 

DIRECTION 

330 

c 

TRIOIR 

OF  the  spherical  triangle  In  THE  FIRST 

GUESS 

PRO- 

335 

c 

CEDURE 

WHENEVER  THE 

ROUTE 

CROSSES  THE  EQUATOR 

• 

336 

c 

337 

h TR  IFG 

=  0 

338 

c 

339 

c 

END  OF 

FUNC  TION  ICONV 

=  1  • 

RETURN  TO  MAIN  PROGRAM. 

340 

c 

42 


RETURN 


************************************* 

*  PART  tl  * 

*  FUNCTION  ICONV  *2  * 

*  eo  crossing  pt  j  angle  * 

************************************* 

determine  WHETHER  point  ‘A*  OR  POINT  *B‘  is  to  serve  as  point 
•£.'  whichever  is  closest  to  the  equator  will  serve  as  *E.‘ 

THE  OTHER  WILL  SERVE  As  *F.‘  NS  I GN  Is  A  FLAG  INDICATING  THE 
ROUTE  INVOLVES  AN  EQUATOR  CROSSING.  NORMALLY  1.  NS  I GN  TAKES 
ON  THE  VALUE  2  FOR  ROUTES  CROSSING  THE  EQUATOR.  NCASE  IS 
A  FLAG  THAT  HAS  THE  VALUE  12  rfHEN  POINT  'A*  SERVES  AS  POINT 
' E •  AND  THE  VALUE  34  WHEN  POINT  'B'  SERVES  AS  POINT  *E.‘ 

NSIGN  =  1 


359 

IF  (FRMLAT 

.LT. 

360 

IF  (FRMLAT 

.GT. 

361 

XL  ATE 

= 

FRMLAT 

362 

XLNGE 

= 

FRMLNG 

363 

PH  I  E 

X 

PHI  A 

364 

XLAtF 

= 

tolat 

365 

XLNGF 

= 

tolng 

366 

PHIF 

= 

PHIB 

367 

NCASE 

= 

1  2 

368 

IF  ( 

ABS(XLATE) 

369 

XLATE 

= 

TOLAT 

370 

XLNGE 

= 

TOLNG 

371 

PHIE 

= 

PHIB 

372 

XLAtF 

= 

FRMLAT 

373 

XLNGF 

= 

FRMLNG 

374 

PHIF 

= 

PHI  A 

375 

NCASE 

34 

205 

C 

C6000 

c 

c 

c 

c 

c 

c 


CONTINUE 

IF  (NOEBUG  .GT.  I>  WRITE  ( 6 • 6000  >  NSIGN,  NCASE 
FORMAT  ( IH  ,  T4.  ‘NSIGN, NCASE  =  ',  T36,  2110) 

ESTABLISH  initial  heading  OF  AIRCRAFT  AT  POINT  ‘E’  on  the  great 
circle,  convert  HEADING  AT  ‘E*  IN  DEGREES  TO  angle  ‘E‘  in 
RADIANS.  TAKE  ITS  TANGENT  TANE  FOR  LATER  USE  IN  THE  FIRST 
GUESS  PROCEDURE  *2. 

HOGE  *  HOG  ( XL ATE »  XLNGE ,  XlATF,  XLNGF) 

IF  (HDGE  .GE.  270.0)  GO  TO  212 
IF  (HOGE  .GE.  1B0.0)  GO  TO  214 
IF  (HDGE  .GE.  90.0)  GO  TO  216 
GO  TO  218 


390 

212 

E  *  HDGE  - 

270.0 

391 

GO  TO  220 

392 

2  14 

E  =  270.0  - 

HDGE 

393 

GO  TO  220 

394 

216 

E  =  HDGE  - 

id 

o 

• 

o 

395 

GO  TO  220 

39t 

2  18 

E  ■  90.0  - 

hdgf 

397 

GO  TO  220 

398 

220 

TANE  =  TAN 

(E  *  (Pl/180.0) ) 

399 

C 

IF  (NDEBUG 

.GT.  1)  WRITE  (6. 

6010) 

HDGE,  E 

400 

C6010 

FORMAT  ( 1H 

,  T  4 ,  ‘HEADING  AT 

POINT 

E  =  '  ,  T  36 , 

401 

c 

•  DEG  * / 1 X , 

T4 ,  T  4 ,  'ANGLE  F 

*  •  t 

T36.  F10.4, 

SET  LATITUDE  TO  ZERO  AND  DETERMINE  LONGITUDF  WHFRE  GREAT  CIRCLE 
CROSSES  EQUATOR.  THIS  IS  POINT  A-PRIME.  WHOSE  LONGITUDE  IS 
XLNGAP  IN  DEGREES. 

XLAT  =  0.0 
GO  TO  500 


230 

C 

C6020 

C 

c 

c 

c 

c 

c 

c 


XLNGAP  =  XLNG 

IF  (NDE8UG  .GT.  1)  WRITE  (6.6020)  XLNGAP 

FORMAT  (1M  .  T  4 .  'LONGITUDE  OF  EQUATOR  CROSSING  =  T36.  F10.4, 

•  OEG • ) 

IF  LONGITUDE  XLNGAP  OF  EQUATOR  CROSSING  POINT  'A-PRIME*  IS 
WEST  OF  LONGITUDE  OF  POINT  'A*  ( FRMLNG ) •  THE  SPHERICAL 

TRIANGLE  OPENS  EASTWARD  (TRIDIR  =  -1.0).  OTHERWISE.  IT 
OPENS  WESTWARD  (TRIDIR  =  1.0). 

XLGAPP  *  XLNGAP 
FRMLNN  «  FRMLNG 

IF  (XLGAPP  .LT.  -90.0)  XLGAPP  =  XLGAPP  ♦  360.0 
IF  (FRMLNN  .LT.  -90.0)  FRMLNN  =  FRMLNN  ♦  360.0 


245 

250 

C 

C6030 

C 

c 

c 

c 

c 

c 

260 


IF  (FRMLNN  -  XLGAPP  )  240.245.245 

TRIDIR  =  -1.0 

GO  TO  250 

TRIDIR  =  1.0 

CONTINUE 

IF  (NDEBJG  .GT.  1)  WRITE  (6,6030)  TRIDIR 
FORMAT  (1H  ,  T  4 ,  'TRIDIR  =  T36,  F10.4) 

establish  initial  heading  of  the  aircraft  at  point  a-prime  on 

CIRCLE.  CONVERT  HEADING  AT  A-PRIME  TO  ANGLE  A-PRIME.  CONVERT 

angle  a-prime  to  radians  and  take  its  sine  for  later  use  in 
first  guess  PROCEDURE  W 3 « 

HDGAP  =  HOG  ( XL  A  T •  XLNGAP,  TOLAT.  TOLNG) 


437 

IF  (HDGAP  .GE. 

270.0)  GO 

TO 

270 

4  38 

IF  (HDGAP  .GE. 

1  BO  *  0  )  GO 

TO 

278 

439 

IF  (HDGAP  .GE. 

90.0)  GO 

TO 

2  80 

440 

GO  TO  285 

44  1 

270 

APRIMD  =  HDGAP 

-  270.0 

442 

GO  TO  290 

443 

275 

APRImD  =  270.0 

-  HDGAP 

444 

GO  TO  290 

4  4  5 

280 

APRIMD  =  HDGAP 

-  90.0 

446 

GO  TO  290 

447 

285 

APRIMD  =  90.0 

-  HDGAP 

448 

GU  TO  290 

290 

C 

C  6  04  0 

c 

c 

c 


SAPRIM  =  SIN  (APRIMD  *  (PI/180.0)) 

IF  (NDEBUG  .GT.  1)  WRITE  (6. 6040)  APRIMD 

FORMAT  (1H  ,  T  4 ,  'APRIMD  =  ',  T36,  F10.4.  •  DEG') 

END  QF  FUNCTION  ICONV  =  2.  RETURN  TO  MAIN  PROGRAM. 


PART  III 

FUNCTION  ICONV  -  3 
XLNG  TO  XLAT 


CONVERT  INTERMEDIATE  LONGITUDE  XLNG  TO  INTERMEDIATE  PHI. 
ICONVX  =  1 

CALL  SPHGLO  (ICONVX,  0.0.  XLNG.  THETAX.  PHI  I 


PART  IV 

FUNCTION  ICONV  =  4 
PHI  TO  XLAT 


CALCULATE  THETA  FROM  PHI. 


c 


477 

478 

479 

480 

481 

482 

483 

484 

485 

486 

487 

488 

489 

490 

491 

492 

493 

494 

495 

496 

497 

498 

499 

500 

501 

502 

503 

504 
506 

506 

507 
50a 

509 

510 

511 

512 

513 
5  1  4 
5  1  6 

516 

517 

518 

519 

520 
52  1 
522 
52  3 

524 

525 

526 

527 
526 

529 

530 

531 
53? 

533 

534 
5  35 

536 

537 

538 

539 

540 
54  1 

542 

543 

544 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


400  THETA  *  ATAN(-0TZ  /  (OTX  »  COSCPHI)  ♦  OTY  *  SIN(PHI)  )  ) 

IF  (THETA!  410,900.900 
410  THETA  =  THETA  ♦  PI 
GO  TO  900 


************************************* 

*  PART  V  * 

*  FUNCTION  ICONV  =5  * 

*  xlat  to  xlng  * 

************************************* 

SAVE  AZIMUTH  ANGLE  OF  THE  CURRENT  POINT  IN  PHIC  FOR  LATER  USE  IN 
CHECKING  FOR  MONOTONIC  PROGRESS  IN  LONGITUDE. 

500  PHIC  =  PHI 

CONVERT  INTERMEDIATE  LATITUDE  XLAT  TO  INTERMEDIATE  THETa. 

I  CON VX  =  1 

CALL  SPHGLO  (ICONVX,  XLAT.  0.0.  THETA,  PMIX) 


************************************* 

*  part  vi  * 

*  FUNCTION  ICONV  =6  * 

*  THETA  TO  XLNG  * 


CALCULATE  PHI  FROM  THETA  (NEKTON’S  ITERATIVE  TECHNIQUE  KITH  AN 
APPROPRIATE  FIRST  GUESS). 

THREE  METHODS  CAN  BE  USED  AS  A  FIRST  GUESS  FOR  THE  AZIMUTH 
ANGLE  phi.  THE  FIRST  And  least  SOPHISTICATED  METHOD  IS 
FOR  THE  USER  TO  SPECIFY  A  FIRST  GUESS  BY  SETTING  INPUT 
IGUESS  =  I  AND  BY  INSERTING  A  FIRST  GUESS  LONGITUDE  IN 
XLNGG  IN  HIS  CALL  TO  THE  GRTCIfl  SUBROUTINE.  THIS  METHOD 
IS  IMPLEMENTED  IN  S  ATEMENTS  605-609. 

THE  SECOND  FIRST  GUESS  METHOD  IS  THAT  OF  APPROXIMATELY 

planar  triangles,  this  method  has  important  limi¬ 
tations  ON  ITS  reliability  and  so  is  used  only  to  find  the 

EQUATORIAL  LONGITUDE  OF  THE  GREAT  CIRCLE  (FUNCTION  ICONV  = 

2).  THE  APPROXIMATELY  PLANAR  TRIANGLE  METHOD 
is  implemented  in  statements  6io-624. 

THE  THIRD  and  MOST  SOPHISTICATE  first  GUFSS  METHOD  IS  THAT 
OF  EXACT  SPHERICAL  TRIANGLES.  HIS  METHOD  REOUIRES 
KNOWLEDGE  OF  THE  EQUATORIAL  C  ING  POINT  AND  ANGLE. 

THE  EXACT  SPHERICAL  TRIANGlE  i  T  GUESAS  METHOD  IS 

implemented  In  statements  625-700. 

625-700. 

600  IF  (IGUESS  .60.  1)  GO  TO  605 
IF  (ICONV  .EO.  2)  GO  TO  610 
GO  TO  625 


*****f irst  guess  method  m***** 

SPECIFIED  FIRST  GUESS.  STATEMENTS  605-609. 

605  ICONVX  *  1 

OUMLAT  =  0.0 

dumthe  =  0.0 

CALL  SPHGLO  (ICONVX,  DUMLAT.  XLNGG.  DUMTHE .  PHIG) 
PHIN  =  PHIG 


84  5  ICONVX  *  2 

546  CALL  SPHGLO  (ICONVX,  DUML  AT ,  XLNGG.  DUMTHE,  PHIG> 

547  609  GO  TO  750 

548  C 

549  C 

550  C  *****FIRST  GUESS  METHOD  #2*4*** 

55\  C 

552  C  APPROX IMATELV  planar  triangle  FIRST  GUESS.  STATE- 

553  C  MENTS  610-624. 

554  C 

555  610  IF  (NSIGN  .Eo.  1)  GO  TO  612 

556  IF  (NSIGN  .Eo.  2)  GO  TO  616 

557  612  PHIN  =  PHIE  +  ABS(  XL  ATE  *  (PI/180.0)  *  TANE  ) 

55b  PHIG  =  PHIN 

559  GO  TO  620 

560  616  PHIN  =  PHIE  ♦  ASS (  XL ATE  *  (PI/160. 0)  /  TANE  ) 

561  PHIG*  PH  IN 

562  GO  TO  620 

563  620  ICONVX  =  2 

564  CALL  SPHGLO  (ICONVX.  OUMLAT,  XLNGG,  DUMTHE.  PHIG) 

565  C  IF  (NDEBUG  .GT.  1)  WRITE  (6,6050)  PHIE,  XLATE.  PHIN 

566  C  605  0  format  (1H  ,  T4.  ' PH  I E . XL  A TE . PH  I N  =  ••  T36.  3F10.4. 

567  C  •  RAO. DEG. RAD1 ) 

568  C  IF  (NDEBUG  .GT.  1)  WRITE  (6.6050)  XLNGG.  PHIG 

569  C6060  FORMAT  (1H  .  T4 .  'FIRST  GUESS  LONGITUDE  A  AZIMUTH  = 

570  C  T36,  2F10.4,  •  DEG. RAO') 

571  624  GO  TO  750 

57?  C 

573  C 

574  C  4****FIRST  GUESS  Ml THOO  «3»***» 

575  C 

576  c  EXACT  SPHERICAL  TRIANGLE  FIRST  GUESS.  STATEMENTS  625-700. 

577  C 

576  C  LATITUDE  XlAT  OF  THE  INTERMEDIATE  POINT  GIVES  THE  ANGULAR 

579  C  DISTANCE  OE-PRlME  IN  RADIANS. 

580  C 

Sol  625  DEPRIM  =  XlAT  *  (PI/180.0) 

58?  C  IF  <  NDE8UG  .GT.  1)  WRITE  (6,6070'  OEPRlM 

583  C6070  FORMAT  (1H  .  T4,  ■ DEPR I M  =  •,  T36,  F10.4.  •  RAD") 

504  C 

565  C  LAW  OF  SINES  FOR  SPHERICAL  TRIANGLES  GIVES  SOLUTION  FOR  ANGULAR 

586  C  DISTANCE  APRMD(l)  AND  APRMD ( 2  )  IN  RADIANS.  APRMD(?>  IS 

587  C  ( I80.0-APRM0( 1 ) )  AND  HAS  ThE  SAME  SINE  AS  APRMO(l).  This  PRO- 

588  C  CEOUPE  IS  NECESSARY  BECAUSE  T HE  ARC  SINE  FUNCTION  RETURNS 

569  C  ONLY  THE  PRINCIPAL  ANGLE. 

590  C 

591  SR  A  T 10  =  SIN(DEPRIM)  /  SAPR1M 

592  C 

593  C  IF  THE  ABSOLUTE  VALUE  OF  S  R  A  T I  .•  EXCEEDS  THE  LEGAL  RANGE.  THEN 

694  C  THE  SOLUTION  IN  LONGITUDE  IS  INDETERMINATE  FOR  A  GIVEN  L  AT  I - 

595  C  TUDE.  IN  SUCH  A  CASE,  THE  METHOD  MAY  6E  SAID  TO  HAVE  FAILED 

596  C  TO  CONVERGE,  AND  ERROR  (ICPNV  *  11)  IS  ASSIGNED. 

5  9  7  C 

59b  IF  I  ABS(SRATIO)  .LE.  1.0)  GO  TO  626 

599  I C  ON  V  =  1  1 

600  GO  TO  1010 

601  626  APRMQ(l)  r  ARSIN  (  SRATIO  ) 

662  APRMJ(?)  =  PI  -  APRMO(l) 

603  C  IF  (nPEBUG  .GT.  1)  WRITE  (6.60-0)  APRMD(l),  APRMD ( 2 ) 

604  C6080  FORMAT  (1H  ,  T4,  • APRMO (  1  )  . APRMD ( 2  )  *  •,  T36,  2F10.4.  '  RAD') 

605  C 

606  c  inverse  of  formula  for  great  circle  distance  (see  subroutine 

pc?  c  'DIstan*  >  gives  first  guesses  xlngggii),  i  *  1.2,  for  latitude: 

60b  C  OF  INTERMTOIATE  POINT  IN  DEGREES.  ThE  FACT  THAT  THERE  ARE 

609  C  T«n  XLNGGC.  VALUES  EMERGES  FRO  i  THE  FACT  THAT  THERE  ARE  TWO 

610  C  APRMD  VALUES.  SIGN  MULTIPLIER  TRIDIR  (-1.0  OR  1.0)  CAUSES 

611  C  ADDITION  OF  LONGITUOE  ELEMENT  FOR  THE  CASE  WHERE  P,INT 

612  C  'A-PRIME'  IS  EAST  OF  POINT  'A'  (SPHERICAL  TRIANGLE  OF-NS 


<4  h 


6  1  3 

614 

615 
6  1  6 
617 
6  1  8 

619 

620 
62  1 
622 

623 

624 

625 

626 
627 
626 

629 

630 

631 

632 

633 

634 
o  35 

636 

637 
t  Jb 

639 

640 
04  1 
042 
643 
o  4  4 
046 

646 

647 
6  4b 
64V 
550 
66  1 

652 

653 
6  5  4 
o65 
o  6  6 
657 
65b 
659 
6  6  0 
66  1 
6  6  ? 
60  3 
0  6  4 
6  O  '■ 
066 
o  6  7 
0  6  H 
069 
6  70 
67  1 
o  7  ? 
6  7  3 
6  7  4 
o  7  5 
6  7  6 

677 

678 
6  79 
680 


C  WESTWARD).  OTHER*  ISE ■  LONGITUDE  ELEMENT  IS  SUBTRACTED. 

C 

C  SPECIAL  CODE  IS  NEEDED  FOR  THE  CASE  OF  AN  EQUATOR  CROSSING. 

C  THE  CODE  1MMC0IATELY  FOLLOWING  AND  EXTENDING  THROUGH  STATE- 

C  MENT  629  ACTS  TO  REVERSE  THE  DIRECTION  OF  THF  SPHERICAL 

C  TRIANGLE  AT  THE  POINT  WHERE  As  AIRCRAFT  PATH  CROSSES  THE 

C  EQUATOR  SOUTHWARD  FROh  THE  NORTHERN  HEMISPHERE  OR  NORTHWARD 

C  FROM  THE  SOUTHERN  HEMISPHERE.  ONLY  IF  nSIGN  EOUALS  2  WILL 

C  AN  EOUATOR  CROSSING  OCCUR.  An  EQUATOR  CROSSING  IS  DECIARED 

C  WHEN  THE  CURRENT  POINT  XLATC  IS  0.0  AND  THE  NEXT-POINT  IS 

C  NEGATIVE  (FOR  A  SOUTHWARD  CROSSING)  OR  POSITIVE  (FOR  A 

C  NORTHWARD  CROSSING).  WHETHER  A  FLIGHT  IS  PROCEEDING  NORTH- 

C  WARD  OR  SOUTHWARD  Is  INDICATED  BY  THfc  XNCRTH  FLAG. 

C 

IF  (NTRIFG  ,E0.  1)  GO  TO  629 

IF  (NSIGN  • E  0 .  1)  GO  TO  629 

IF  (XNORTH)  627,628,628 

627  IF  (XLATC  .NE.  0.0)  GO  TO  629 

IF  ( XL AT  .GE.  0.0)  GO  TO  629 

TRIDIR  =  -TRIOIR 

NTRIFG  =  1 

GO  TO  629 

628  IF  (XLATC  .NE.  0.0)  GO  TO  629 

IF  ( XL  AT  .LE.  0.0)  GO  TC  629 

TRIDIR  =  -TRIOIR 

NTRIFG  =  I 

629  CONTINUE 

C  IF  (NOEHUG  .GT.  1)  WRITE  (6.6090)  TRIDIR 

C6090  FORMAT  (1H  .  T4,  'TRIDIR  =  T36,  F10.4) 

C 

DO  640  I  =  1.2 

630  XLnGGG(I)  =  XLNGAP  ♦  TRIDIR  ♦  (  (180.0/PI)  ♦ 

X  ARCOS  (  COS ( APRMD( I ) )  /  COS(  * L A T  *  (PI/180.0)  )  )  > 

c 

C  ADJUST  COMPUTFD  LONGITUDE  XlNGo  FOR  DATELINE  FOLD. 

C 

IF  (XLNGGG(I)  .GT.  180.0)  GO  TO  632 
IF  (XLNGGG(I)  .LT.  -180.0)  GO  TO  634 
GO  TO  636 

632  XLNGGG(I)  =  XLNGGG(I)  -  360.0 
GO  TO  636 

634  XLNGGG(I)  =  XLNGGG(I)  ♦  360.0 
GO  TO  636 
C 

C  CONVERT  LONGITUDES  XLNGGG(I)  IN  DEGREES  TO  AZIMUTHS  PHIGG(I)  IN 

C  RADIANS. 

C 

036  ICONVX  *  1 

CALL  SPHGLO  (ICONVX,  DU“LAT,  XLnGGG(I),  DUMTHF.  PHIGG(I)) 

C  IF  (NDEHUG  .GT.  1)  WRITE  (6,6100)  XLNGGG(I) 

CftlOu  FORMAT  (1H  .  T  4  .  'XLNGGC-(I)  =  *,  T36.  Flo. 4,  •  DEG') 

640  CONTINUE 

C 

C  INITIALLY,  WHICHEVER  LONGITUDE  XLNGGG(I)  HAS  AN  ASSOCIATED 

C  AZIMUTH  PHIGG(I)  ABSOLUTELY  CLOSER  TO  PHIC,  THF  AZIMUTH 

C  ANGLE  OF  THE  CURRENT  POINT,  13  SELECTED  AS  THE  FIRST  GUESS 

C  LONGITUDE  XLNGG  IN  DEGREES  WITH  CORRESPONDING  AZIMUTH  PHIG 

C  (POINTER  IS  SAVED  IN  IGOOD).  IF  XE  AST  IS  GREATER  THAN 

C  ZERO  (Ltss  THAN  ZERO).  PHI^,  MOST  BE  GREATER  THAN  UESS 

c  THAN)  PHIC.  IF  THIS  is  NOT  The  CASE.  THE  ALTERNATE  L0N61- 

C  TUDE  SOLUTION  IS  USED  AS  FIRST  GUESS  (POINTERS  IGOOD  AnD 

C  1 8  AD  ARE  SWAPPFO). 

C 

C  IF  (NOE0UG  .GT.  1)  WRITE  (6,6110)  PHIGG.  PuIC.  XEAST 

Co  1  1  0  format  (  1 H  ,  T4,  ' PHIG (  1  ) -AND- ( 2 > .PHIC , XEAST  *  ',  T36,  4F10.4) 
IGOOD  =  1 
I  BAD  =  2 

IF  (  ( ABS(PHIGG(  IBAD)  -  PHIC))  .GT.  ( AB S ( PH  I GG ( I  GOOD )  -  PHIC))  ) 


■iC i, -  - 


Ob  1 
6  82 
o  b  3 
0B4 
685 
066 
06  7 
o  b  8 
o89 
to  9  0 
o  9  1 

6  9? 

693 

694 
69*. 

696 

697 

698 

699 
70  C 
701 
70? 
703 

70  A 
70S 

7  u6 

707 

708 

709 
7  1  C 
7  1  1 
71? 

71  3 
7  1  4 
7  1  5 
V  1  6 
71  7 
7  l  8 
7  1  9 
7  20 
71  1 

let 
7  2  3 
724 
7r  S 
7  2  6 

72 / 

7  2  8 
729 
MO 

73  1 
7  3c 
733 
7  34 
7  3  0 
7  3  6 
737 
736 

739 

740 

74  1 
74? 


& 


645 

650 

652 

654 

660 

C 

C 

C  6  1  20 
C 

c 

c 

c 

c 

c 

c 

690 

700 

C 

C 

C 

C 

c 

c 

c 

c 

750 


C 

c 

c 

775 

C 

c 

c 

c 

0 

c 

v. 

c 

c 

t 


800 


850 


-  PHIC)  650,652*652 
660,654 , 654 
654,654,660 


GO  TO  645 
ITEMP  *  IGOHO 
IGOOD  =  I B  A  0 
I  3  A  D  *  ITEMP 
IF  ( PHlGGt IGOOD ) 

IF  (XEAST) 

IF  (XEAST) 

ITEMP  *  IGOOD 
IGOOD  *  IBAD 
IBAQ  =  ItEMP 
CONTINUE 

IF  (NDEPUG  .GT.  1)  WRITF  (6,61?0) 

XLNGGG (  I  BAD > ,  PHIGG(IBAD) 

format  (1H  ,  T  4 ,  * XLNGGG, PHIGG  FOR  IGOOD  =  •,  T36, 
T4.  * XLNGGG, PHIGG  FOR  IBAD  =  »,  T36,  2F10.4) 


XLNGGG (  I  GOOD  )  ,  PH  I GG (  I  GOOD ) , 


2F10.4/1X, 


SET  FIRST  GUESS  LONGITUDE  AnD  PHI. 


XLNGG  *  XLNGGG  (  I  GOOD ) 
PHIG  =  PHIGG  (IGOOD) 

STORE  PHIG  IN  PHIN. 


PHlN  =  PHIG 
GO  TO  750 


♦****Nt WTON • S  ITERATIVE  TECHNIQUE***** 

imppjve  by  newton's  method. 

INITIALIZE  iteration  count  and  DENOMINATOR. 

N  =  0 

NP 1  =  N  ♦  1 

TANTMF  =  TAN(ThETA) 

F  UNCT I  ON . 

FPhI  =  GTX  *  CCS(PHIN)  ♦  QTY  *  SIN(PHIN)  ♦  Q  T  2  /  TANTHE 

FIRST  DERIVATIVE  OF  FUNCTION. 

Fprnl  =  - 0  T  X  *  SIN (PHIN)  ♦  QTY  *  COS( PH  IN) 

Improved  phi. 

PhlNPl  =  PHlN  -  ( F  P  H 1  /  FPOHI  ) 

TEST  FOR  CONVERGENCE.  STOP  WITH  ERROR  CONDITION  ICDNV  =  11  IF 

METHOD  FAILS  TO  CONVERGE  IN  N  H  A  X  ITERATIONS. 

DELPHI  =  AdStPHlNPl  .  PHIN) 

IF  (OFLPHI  -  EPS)  850*850, HOO 
PHlN  =  PH  I NP 1 
N  =  N  ♦  1 

NP  1  =  N  ♦  1 

IF  (N  .L  T .  h  M  A  X  )  GO  TO  775 
ICONV  =  ll 
GO  TO  1010 
PHI  =  PH  I  NP 1 


C 


743 

744 

745 

746 

747 
74P 

749 

750 

751 

752 

753 

754 

755 

756 

757 

758 

759 

760 

761 

762 

763 

764 

765 

76  6 
767 
760 
769 
7  70 

77  1 

772 

773 

774 

775 

776 

777 

778 

779 
700 
7b  i 
782 
7  M  3 
784 
70  5 
786 
767 
788 


C 

C 

c 

c 

c 

c 

c 

c 

900 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

1000 

6  1  30 

& 

& 


1010 

1020 

it 

1050 
1  060 

1070 


*  PART  VII 

*  ALL  FUNCTIONS 

*  CONVERT  SOLUTION  (THETA, PHI) 

*  sack  to  global  <xlat,xlng> 

♦**♦♦***♦*♦♦♦♦**♦♦♦♦♦*♦♦♦♦♦♦*♦***♦♦♦ 


I  CON VX  =  2 

CALL  SPHGLO  (ICONVX*  XLAT.  XLNg*  ThETA.  PHI) 
IF  (  I  CON V  .EG.  2)  GO  TO  230 


PART  VIII 

NORMAL  termination 


RETURN 


PART  IX 

ABNORMAL  termination 


WRITE  (6.6130)  FRMLAT,  FRMLNG.  TOLAT,  TOLNG 

format  <iho.  •♦♦♦♦wind  factor  r  r  r  o  r  di  agnostics^^  »/i  x , 

•FROM-LAT  =  »,  F 1 0 • 2 ♦  T30.  • F  R  OM- L  NG  s  ••  T41,  F10.2/1X. 

• TG-L  A  T  =  •.  F 1 0 • 2 .  T32.  *  TO-LNG  =  T41.  F10.2//) 

CALL  FXEM  <61.  'SUBROUTINE  GRTC I R  .  •  •  I LLE G Al  ICONV,  9> 

GO  TO  1050 
CONTINUE 
GO  TO  1050 

WRITE  (6.6130)  FRMLAT.  FRMLnG.  TOLAT,  TOLnG 

CALL  FXEM  (61.  ‘SUBROUTINE  GRTC I R . . . SOLUT  I  ON  NOT  UNIQUE  IN  THETA 
12) 

IF  (ICONVS  • E  0 •  2  .OR.  ICONVS  .EQ.  3)  GO  TO  1060 

IF  (ICONVS  .EQ.  4  .OR.  ICONvS  .EQ.  5)  GO  TO  1070 

XL  A  T  =  0.0 

THETA  =  0.0 

RETURN 

XLNG  =  0.0 

PHI  =  0.0 

return 

END 
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1  CSPMGLO  CON V  SPH<> GL OB AL /R  •  C.  WHIT0-./05  FEB  1979 

2  C 

3  SUBROUTINE  SPHGlCj  (ICONV,  XL*T,  XLON.  THETA,  PHI) 

4  C 


6 

C* 

7 

C* 

program 

ID- 

SPHGLO 

8 

O 

met  analyst- 

MAJ  ROGER 

C.  wHITON.  USAFETAC/DNS.  EXT  5*12 

9 

c* 

sys  analyst- 

MAJ  ROGER 

C.  WHITON.  USAFETAC/DNS.  EXT  5*12 

10 

c* 

programmer- 

MAJ  ROGER 

C.  WHITON.  USAFETAC/DNS.  EXT  5*12 

1  1 

c* 

1  2 

c* 

CREATED 

1 

2 

o 

OS  FEB  197S  PROJECT-  192301 

1  3 

c* 

14 

c* 

DESCRIPTION- 

this  sooroutime  su3prog»am  convfrts  global  coor¬ 

1  5 

c* 

dinates  <LAT.t.0N)  To  MATHEMATICAL  SPHERICAL  COOR¬ 

1  6 

c* 

DINATES  (THETA, PHI ),  OR  VICE-VERSA,  DEPENDING  ON 

17 

c* 

THE  VALUE 

OF  THE  ICONV  FLAG  (ICONVsl  FOR  CONVER- 

1  ft 

c* 

SION  OF  (LAT.LON)  TO  (THETA, PHI).  AND  IC0NV*2  FOR 

19 

c* 

CONVERSION 

OF  (TrifTA.PHI)  TO  (LAT.LON). 

20 

c* 

21 

c* 

ME  THO0- 

FOR  CONVERSION  OF  (LAT.LON)  TO  ( T HE T A , P M I ) . . . 

22 

c* 

23 

c* 

THETA  * 

( P I /2 . 0  )  -  XL  AT 

24 

c  * 

PHI 

-XLON  FUfc  XLON  ,LE.  0.0 

25 

c* 

PHI  = 

<2.0*EI)  -  XLON  FOR  XLON  .GT.  0.0 

2 6 

c* 

IF  (PHI 

.GT.  PI)  PHI  a  PHI  -  ( 2  •  0  *  P 1  ) 

27 

c* 

26 

c* 

FOR  CONVERSION  Or  (THETA, PHI)  TO  (LAT.LON)... 

29 

c  * 

30 

c* 

XL  AT  a 

(PI/2.0)  -  THETA 

3  1 

c* 

XLON  = 

-PHI  FOR  PHI  • L  T .  PI 

32 

c* 

XLON  = 

( 2 • 0*P I  )  -  PHI  FOR  PHI  .GE.  PI 

33 

o 

34 

c* 

FOR 

ABS(PHI)  .LT.  2  »  0  *  P I 

35 

c* 

36 

c  * 

NOTE  THAT 

SOUTH  LATITUDES  and  EAST  LONGITUDES  ARE 

37 

c* 

HANDLED  AS 

NEGATIVE  NUMBER'  AND  MUST  BE  SO  SPECI- 

36 

c* 

FIED  ON  INPUT/OUTPUT  FROM  THE  SUBROUTINE.  L  A  T I - 

39 

c* 

TUOES  AND 

LONGITUDES  ARE  IN  DEGREES,  WHILE  THETA 

40 

c* 

AND  PHI  ARE  In  RaDIAnS.  Thf  AZIMUTH  ANGLF  PHI  IS 

4  1 

c* 

THE  ANGLE 

COUNTERCLOCKWISE  FROM  GREFNWICH.  THE 

4  c 

C* 

ABSOLUTE  VALUf  BINGES  FROM  0.0  TC  E.O*Pl  RADIANS. 

43 

c* 

NEGATIVE  VALUf S  INDICATE  CLOCKWISE.  THF  COLATI- 

44 

c  * 

TUDE  ANGLE 

THETA  is  the  ANGLE  SOUTHWARD  FROM  THE 

45 

c* 

NORTH  POLE 

and  ranges  from  o.o  to  pi  radians. 

46 

c* 

47 

c* 

input  n 

ICONV  = 

FLAG  CONTROLLING  WHETHER  CONVERSION  IS 

4  8 

c  * 

OUTPUT- 

FROM  (LAT.LON)  TO  (THETA, PHI)  (1C0NV=1> 

49 

c* 

OR  FROM  (THETA, PHI)  TO  (LAT.LON) 

SC¬ 

c* 

(ICl)NVs?).  WHEN  RFTURNFD  AS  OUTPUT, 

SI 

C  4 

ICONV  -  1  OR  ?  FOR  NORMAL  TFRM] NA¬ 

52 

c  * 

TION  0'«  A  NUMBER  GRfrATLR  THAN  OR 

5  3 

c  * 

OR  r  Q U  \ L  TO  ID  FOR  ABNOPMAL  TERMINA¬ 

54 

c* 

TION.  ABNORMAL  TERMINATION  codes 

5  5 

c* 

FOR  ICuNV  ARE •  •  • 

S  6 

c* 

57 

c  * 

10-  ERRONEOUS  ICONV  SPECIFIED 

58 

c  ♦ 

11-  OUT-OF-BOUNDS  XLAT 

39 

C4 

12-  OUT -OF -HOUNDS  XLON 

60 

c  ♦ 

13-  OUT-OF-BOUNDS  THETA 

6  1 

C* 

14-  OUT-OF-BOUNDS  PHI 

62 

c* 

6  3 

C  ♦ 

XL  AT  = 

LATITUDE  IN  DEGREES,  0.0  TO  90.0, 

54 

c* 

PQSITIVF  FOR  NORTH  LATITUDES  AND  NEGA¬ 

65 

C* 

TIVE  FUR  SOUTH  LATITUDES. 

66 

c  * 

67 

c* 

XLON  = 

LONGITUDE  IN  DEGREES.  0.0  TO  IbO.J,  POSI¬ 

68 

c* 

TIVE  Flu  WtST  LONGITUDES  AND  NEGATIVE  FOR 

EAST  LONGITUDES 


69 

C* 

70 

C* 

71 

c* 

72 

c* 

73 

c* 

74 

c* 

75 

c* 

76 

c* 

77 

c* 

78 

c* 

79 

c* 

system  sub¬ 

80 

c* 

programs 

ai 

c* 

USED- 

82 

c  * 

83 

c* 

USER  SUB¬ 

34 

c* 

PROGRAMS 

b5 

c* 

USED- 

66 

c* 

87 

c* 

ESTImATEO 

88 

c* 

CPU  TIME- 

69 

c* 

90 

c* 

STORAGE 

9  1 

c* 

REOUIREmEnTS- 

92 

c* 

93 

c* 

94 

c  * 

PROGR Am 

95 

c* 

UPDATES- 

96 

o 

THETA  =  COLATITUDE  ANGLE  IN  RADIANS.  0.0  TO  PI 
RADIANS  SOUTHWARD  FROM  NORTH,  POSITIVE. 

PHI  =  AZIMUTH  ANGLE  IN  RADIANS,  0.0  TO  PI 

RADIANS,  COUNTERCLOCKWISE  FROM  GREENWICH 
LOOKING  DOWN  AT  THE  NORTH  POLE.  NEGATIVE 
VALUES  INDICATE  CLOCKWISE. 


NONE 


NONE 


TO  BE  SUPPLIED. 


procedure  plus  data  occupy  isa  words  of 

CORE  STORAGE. 


NONE 


9rt  C 

99  DATA  TWGPI/6. 2831853/,  P I / 3 . 1 4 1 5 9 2 7/ ,  PID2/1. 5707963/, 

100  4  DTOR/O.O 1 745329/ ,  K  T  (,D/ 5  7 . 295  780/ 

101  C 

102  C  DEPENDING  ON  THE  VALUE  OF  THE  ICDNV  FLAG,  CONVERT  (LAT.LON) 

103  C  To  (THETA. PHI)  <IC0NV=1>  OR  CONVERT  (THETA, PHI)  To 

104  C  (LAT.LON)  (IC0NV=2).  FOR  INPUT  1C0NV  OTHER  THAN  l  OR  2, 

105  C  AN  ERROR  CONDITION  (ICONv=10)  IS  RETURNED. 

106  C 

107  IF  (  I C  ON  V  •  F  0 .  1  .OR.  ICONV  .  E •  2)  GO  TO  100 

10H  ICONV  =:  10 

IDS  GO  TO  1000 

110  100  GO  TO  (200,  400),  ICONV 

111  C 

11/  C  CONVERT  (LAT.LON)  TO  (THETA, PHI).  RETURN  FRROR  CONDITION 

I1J  C  IC0NV=11  FOR  OUT-OF-BOUNDS  XLaT  AND  I  CON V= 1 2  FOR  OUT-OF- 

114  C  30UN0S  XLON  . 

ns  c 

116  200  IF  (XLAT  .LF.  90.0  .AND.  XLAT  .GE.  -90.0)  Gn  TO  220 

1 1 7  ICONV  =  1  1 

lib  GO  TO  1000 

119  220  IF  ( X  L  ON  ,LE.  18  0.0  .AND.  XL  ON  .GE.  -180.0)  GO  TO  250 

120  ICONV  =  1 2 

121  GO  TO  1000 

122  250  THETA  =  PID?  -  (XLAT  *  DTOR  ) 

123  IF  (XLDN)  260.260.275 

124  260  PHI  =  -(>  ON  *  DTOR) 

125  IF  (PHI  .GT.  FI)  PHI  =  PHI  -  TWOPI 

126  GO  TO  BOO 

127  275  PHI  =  TWOPI  -  ( XLON  *  DTOR) 

128  IF  (PHI  .GT.  PI)  PHI  =  PHI  -  TwOPl 

129  GO  TO  SOU 

130  C 

131  C  CONVERT  (THETA, PHI)  TU  (LAT.LON).  RETURN  ERROR  CONDITION 

132  C  I C  ONV  =  1 3  FOR  OUT-OF-BOUNDS  THETA  AND  IC0NV=14  FOR  OUT-OF- 

133  C  BOUNDS  PHI.  ADJUST  COMPUTFD  LONGITUDE  FOR  DATELINE  FOLD. 

1  34  C 

1  35 
136 


400  IF  ( AS  S ( THE  T  A )  .LE.  PI)  GO  TO  420 
ICONV  =  13 


137 

GO  TO  1000 

1  38 

420 

IF  (ABS(PHI)  .LE.  T*fOPI)  GO 

1  39 

ICONV  =  14 

14  0 

GO  TO  1000 

141 

450 

XLAT  =  ST  00  *  (PI02  - 

theta  ) 

142 

IF  (PHI  .LT.  PI)  GO  TO 

475 

143 

XLON  =  RTOD  *  (TWOPI  * 

PHI  ) 

144 

GO  TO  500 

145 

475 

XLON  =  -(PHI  *  RTOD) 

146 

500 

IF  (XLON  .GT.  160.0) 

XLON 

147 

IF  (XLON  .LT.  -160.0) 

XLON 

146 

c 

149 

c 

NORMAL  TERMINATION. 

150 

c 

151 

800 

RETURN 

152 

c 

153 

c 

ABNORMAL  TERMINATION. 

154 

c 

155 

1  000 

xlat  =  0.0 

1  56 

XLON  =  0.0 

157 

THEtA  =  0.0 

156 

PHI  =  0.0 

1  59 

RE  TURN 

160 

END 

TO  450 


=  XLON  ♦ 


360 1 
360. 


CD  I  ST  AN  CRT  CIRCLE  DIST/R.  C.  WHITON/26  OCT  1973 


2 

C 

3 

SUBROUTINE 

4 

C 

5 

C+  +  * 

************** 

6 

c* 

7 

c* 

PROGRAM  ID- 

8 

c* 

MET  ANALYST- 

9 

c* 

SYS  ANALYST- 

1  0 

c* 

PROGRAMMER- 

1  1 

c* 

1  2 

c* 

created  on- 

1  3 

c* 

1  4 

c* 

DESCRIPTION- 

1  5 

c* 

1  fc 

c  + 

1  7 

c  ♦ 

1  b 

c  * 

1  9 

c* 

METHOD- 

20 

c* 

2  1 

c* 

22 

c* 

c  J 

c* 

24 

c* 

25 

c* 

26 

c* 

27 

c* 

?b 

c* 

29 

c* 

3  0 

c* 

31 

c* 

32 

c* 

LIMITATIONS- 

3  3 

c« 

34 

c* 

35 

c» 

REFERENCES- 

36 

c* 

37 

c* 

3H 

c» 

39 

c* 

INPUT- 

40 

c« 

4  1 

c» 

4  2 

c* 

4  3 

c* 

44 

c* 

45 

c* 

46 

c* 

47 

c* 

48 

c# 

49 

c* 

OUTPUT- 

50 

c* 

51 

c* 

52 

c* 

SYSTEM  SUD- 

53 

c* 

PROGRAMS 

54 

c* 

USED- 

5  5 

c* 

56 

c  * 

USER  SUB- 

57 

c* 

PROGRAMS 

5  8 

c* 

USED- 

59 

c* 

60 

c* 

ESTIMATED 

6  l 

c* 

CPU  TIME- 

62 

c* 

63 

c* 

64 

c* 

STORAGE 

65 

c* 

REQUIREMENTS 

66 

c* 

DISTAN 

MAJ  ROGER  C.  WHITON,  US AFE T AC /DNS ,  EXT  3412 
MAJ  ROGER  C.  WHITON.  US AF E T AC/DNS ,  FXT  5412 
MAJ  ROGER  C.  WHITON,  US AFE T AC /DNS ,  EXT  5412 


26  OCT  1978 


PROJFCT- 


THIS  SUBROUTINE  SUBPROGRAM  COMPUTES  GREAT  CIRCLE 

distance  over  sphere  of  the  earth  between  two 

POINTS  •  A  •  AND  <B'  WHOSE  (  L  AT  I  TUDE  ,  L  ONG  I  TUDE  )  IS 

specified  in  input  arguments. 

THE  CONVENTIONAL  FORMULA  FOR  THE  GREAT  CIRCLE  ARC 
GCD  IN  RADIANS  IS  USED  (REF  2).  NAMELY. 

GCD(RADIANS)  =  ARCOS ( S IN ( FRMLAT ) *S INC TOL AT )  + 
COS(FRMLAT>*COS<TOLAT)*COS(TOLNG-FRMLNG>) 

where  independent  variables  are  as  described  in 
the  input  arguments  below,  later,  gcd  in  radians 

OF  ARC  IS  CONVERTED  TO  GCD  IN  NAUTICAL  MILES  ( NM ) 

OF  CIRCULAR  DISTANCE  BY  MULTIPLYING  PY  THE 

radius  of  thf  earth,  3440  nm.  gcd  is  returned  to 

THE  MAIN  PROGRAM  IN  NM. 

1.  (LAT.LON)  of  destination  must  not  EE  SAME 
AS  (LAT.LON)  OF  ORIGIN. 

1.  CONRAD  AND  POLLACK,  1950-  METHODS  In  CLIMA¬ 
TOLOGY 

2.  HEWLETT-PACKARD .  1975-  Hp-65  NAVIGATION  PAC  1 

FRMLAT  =  LATITUDE  OF  POINT  'A*  (ORIGIN)  IN  DEG 
FPMLNG  =  LONGITUDE  OF  POINT  •  A  '  (ORIGIN)  IN  DEG 
TOLAT  =  LATITUDE  OF  POINT  'B'  (DESTINATION)  IN 
DEG 

TOLNG  =  LONGITUDE  OF  POINT  'B'  (DESTINATION)  IN 
DEG 

note-  all  south  latitudes  and  east  longitudes  must 
BE  SUPPLIED  as  negative  NUMBERS. 

GCD  =  GREAT  CIRCLE  DISTANCE  IN  NM  OF  ARC  OVER 

THE  SURFACE  OF  THE  EARTH 


SIN,  COS.  ARCOS 


ON  A  OEC-IO  GENERAL  PURPOSE  COMPUTER.  THIS  SUB¬ 
ROUTINE  REQUIRES  2.3  CPU  MILLISECONDS  PER  CALL. 


PROCEDURE  PLUS  DATA  OCCUPY  9A  WOROS  OF 
STORAGE. 


r  ji-  ■  koLv*,*' 


c* 

C*  PROGRAM 

C*  UPDATES-  NONE 

C********************************************************************* 

C  CONSTANTS.  otor  converts  degrees  to  radians,  a  is  the  radius  of 

c  THE  earth  in  nautical  miles  (NM>. 

c 

DATA  DTOR/0. 01745329/.  A/3440.0/ 

c  COMPUTE  GCO  IN  RADIANS  OF  ARC.  MULTIPLICATION  BY  DTOR  CONVERTS 

C  DEGREES  TO  RADIANS. 

C 

SIN1  =  SINCFRMLAT  *  DTOR) 

SIN2  =  S I N  <  T  OL  AT  *  DTOR) 

COS)  =  COSCFRMLAT  *  DTOR) 

COS?  =  COS  C  TOLAT  *  DTOR) 

C0S3  =  COS  <  < TOLNG  *  DTOR)  -  ( FRMLNG  *  DTOR)  ) 

GCR  =  S I N I  *  SIN2  ♦  C0S1  *  C0S2  *  C0S3 
GCO  *  ARCOS(GCR) 

C  CONVERT  GCO  FROM  RADIANS  OF  ARC  TO  CIRCULAR  DISTANCE  BY  MULTIPLY- 

C  ING  BY  THE  RADIUS  OF  THE  EARTH  A,  3440  NM. 

C 

GCO  =  GCO  *  A 
C 

C  TERMINATION. 

C 

RETURN 

END 


CHDG 

C 


INITIAL  HEAOING/P.L.  HEROD/12  FEB  1979 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 
1  1 
12 

13 

14 

15 
If 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 
4  1 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 
5b 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 


C 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

c* 

c* 

c* 

c» 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

C  4 

c* 

c* 

c* 

c* 

c» 

C  4 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c» 

c* 

c« 

c» 


FUNCTION  HDG  < 

****************** 

PROGRAM  ID- 
MET  AnALYST- 
SYS  ANALYST- 
PROGRAMMER- 

CREATED  ON- 

DESCRIPTION- 


METHOO- 


limitations- 


reference- 


INPUT- 


OUTPUT- 


SYSTEM  SUB¬ 
PROGRAMS 
USED- 

USER  SUB¬ 
PROGRAMS 
USEO- 

ESTIMATED 


XL AT  1 •  XLNG1.  XL  A  T2 •  XLNG2 ) 

HOG 

CAPT  PATRICK  L.  HEROD,  US AF ET AC/DnS ,  EXT  5412 
MAJ  ROGER  C.  WHlTON.  USAFET AC/DNS .  EXT  5412 
CAPT  PATRICK  L.  HEROD.  USAFETAC/DNS.  EXT  5412 


12  FEB  1979 


PROJECT- 


192301 


this  function  subprogram  finds  the  initial  heading 
(HDG)  along  a  great  circle  course  flown  from  point 

<A«  ( XL  AT  1 , XLNG 1 )  TO  POINT  'B*  ( XL AT 2 . XLN62 > .  IT 

WAS  DEVELOPED  from  THE  HEWLETT-PACKARD  65  (HP-65) 
AVIATION  PAC  1,  GREAT  CIRCLE  NAVIGATION.  APRIL 
1975. 

THIS  FUNCTION  SUBPROGRAM  FIRST  USES  INPUT  LATITUDES 
AND  LONGITUDES  FOR  POINTS  'A*  AND  *8*  AND  CONVERTS 
THEM  FROM  DEGREES  TO  RADIANS.  IT  THEN  USES  THE 
HEADING  FORMULA  FOUND  ON  PAGE  54  OF  THE  HP-65  AVIA¬ 
TION  PAC  1  FOR  CALCULATION  OF  THE  INITIAL  HEADING 
FROM  POINT  'A.*  SINCE  THE  GREAT  CIRCLE  DISTANCE 
(GCD)  IS  IN  THIS  HEADING  FORMULA.  THE  FUNCTION  SUB¬ 
PROGRAM  MAKES  A  CALL  TO  SUBROUTINE  DISTAN,  WHICH 
RETURNS  THE  GCD  IN  NAUTICAL  MILES  <NM>  FROM  POINT 
•A*  TO  POINT  'B.*  FUNCTION  HDG  CONVERTS  GCD  IN 
NM  TO  GCD  IN  RADIANS.  THE  GCD  IN  RADIANS  IS  USED  IN 
THE  HDG  FORMULA.  HDG  IS  COMPUTED  IN  RADIANS. 

THERE  IS  A  FINAL  CONVERSION  PROCESS  WHICH  ENABLES 
THIS  FUNCTION  TO  OUTPUT  THE  INITIAL  HEADING  IN 
DEGREES. 

TRUNCATION  AND  ROUNDOFF  ERRORS  OCCUR  WHEN  POINT  *A« 
ANO  POINT  'B'  ARE  VERY  CLOSE  TOGETHER  (1  MILE  OR 
LESS).  INPUT  DATA  IS  IN  DECIMAL  DEGREES.  NOT  DE¬ 
GREES.  MINUTES  AND  SECONDS.  NORTH  LATITUDES  ANO 
WEST  LONGITUDES  ARE  POSITIVE  NUMBERS.  SOUTH  LATI¬ 
TUDES  and  east  longitudes  are  negative  numbers. 

HEWLETT-PACKARD,  1975-  AVIATION  PAC  1, 

GREAT  CIRCLE  NAVIGATION 

XLATl  =  LATITUDE  OF  POINT  ' A •  (ORIGIN)  IN  DECIMAL 
DEGREES 

XLNGl  =  LONGITUDE  OF  POINT  »A'  (ORIGIN)  IN  DECIMAL 
DEGREES 

XL AT  2  =  LATITUDE  OF  POINT  • B '  (DESTINATION)  IN 

DECIMAL  DEGREES 

XLNG2  =  LONGITUDE  OF  POINT  »B'  (DESTINATION)  IN 
DECIMAL  DEGREES 

HDG  =  INITIAL  COURSE  HEADING  ALONG  GREAT  CIRCLE 
FROM  POINT  ’A'  TO  POINT  *B*  IN  DECIMAL 
degrees  toward  which 


SIN.  COS.  ARCOS 


DISTAN 


55 


CPU  TIME 


TO  BE  SUPPLIED 


69 

C  4 

CPU  TIME-  TO  BE  SUPPLIED. 

70 

c* 

71 

c* 

STORAGE 

72 

c* 

REQUIREMENTS-  PROCEDURE  AND  DATA  OCCUPY  106  WORDS  OP  CORE 

73 

c* 

STORAGE. 

74 

c* 

75 

c* 

PROGRAM 

76 

c* 

updates-  none 

77 

c* 

78 

c*********«************************************************************ 

79 

c 

80 

OATA  A/3440.0/.  P I / 3 . 1 4 1 592 7/ ,  DT0R/0 . 0 1 74 5329/ 

81 

c 

82 

c 

CONVERT  LATITUDE  AND  LONGITUDE  FROM  DEGREES  TO  RADIANS. 

83 

c 

84 

XLATIR  =  XL AT  1  *  OTOR 

85 

XLNG1R  =  XLNG1  *  DTOR 

86 

XLAT2R  =  XLAT2  *  DTOR 

87 

XLNG2R  s  XLNG2  4  DTOR 

88 

c 

89 

c 

CALCULATE  GCD  IN  NAUTICAL  MILES  (GCDNM)  FROM  C XL AT  1 , XLNG 1 )  TO 

90 

c 

(XLAT2.XLNG2). 

91 

c 

92 

CALL  0 1  ST  AN  (XlATI.  XLNG 1 .  XLAT2.  XLNG2 .  GCDNM) 

93 

c 

94 

c 

CONVERT  GCD  IN  NAUTICAL  MILES  (GCDNM)  TO  GCD  IN  RADIANS  ( GCDR ) . 

95 

c 

96 

GCDR  s  GCDNM  /  A 

97 

c 

98 

c 

CALCULATE  HEADING  IN  RADIANS  (HDGR)  FROM  HEWLETT-PACKARD  65 

99 

c 

100 

HDGR  =  ARCOS< <SIN(XLAT2R)  -  CCOS(GCDR)  *  S I N  1  XL AT l R >  )  )  / 

101 

s 

(SIN(GCOR)  *  COS ( XL AT  1 R  > )  ) 

102 

IF  (SIN(<XLNG1  -  XLNG2 >  *  DTOR)  «LT.  0.0)  HDGR  =  2.0  *  PI  -  HDGR 

103 

c 

104 

c 

CONVERT  HEADING  IN  RADIANS  (HDGR)  TO  HEADING  IN  NAUTICAL  MILES 

105 

c 

( HDG ) . 

106 

c 

107 

HOG  *  HDGR  /  OTOR 

108 

c 

109 

c 

NORMAL  TERMINATION. 

1  1  0 

c 

1  1  1 

RETURN 

1  12 

END 

56 


■1 


1  CBRLNG  BRACKET  LONGITUOE/P.  L •  HEROD/17  JAN  1979 


2 

C 

3 

SUBROUTINE 

SRLNG  (XLNG.  XLNGLO.  XLNGHI) 

4 

c 

5 

c*****««*»**************************************»********»********** 

6 

c* 

7 

c* 

program  id- 

brlng 

8 

c* 

MET  ANALYST- 

CAPT  PATRICK  L.  HEROD.  USAFE TAC/ONS .  EXT  5412 

9 

c* 

SYS  ANALYST- 

MAJ  ROGER  C.  mHITON.  US AFE T AC/DNS .  EXT  5412 

10 

c* 

PROGRAMMER- 

CAPT  PATRICK  L.  HEROD.  US AFE T AC/DNS .  EXT  5412 

1 1 

c* 

1  2 

c* 

CREATED  ON- 

17  JAN  1979  PROJECT  192301 

13 

c* 

14 

c* 

DESCRIPTION- 

THIS  SUBROUTINE  SUBPROGRAM  FINDS  THE  LONGITUDINAL 

1  5 

c* 

VALUES  that  BRACKET  A  GIVEN  LONGITUDE.  THE  LON- 

16 

c* 

GITUOINAL  BRACKET  VALUES  ARE  IN  INCREMENTS  OF  30.0 

1  7 

c* 

I.E..  A  GIVEN  LONGITUDE  OF  132.35  HAS  BRACKET  VALUES 

IB 

c* 

QF  120.0  AND  150.0. 

19 

c* 

20 

c* 

METHOD- 

THE  METHOD  USED  TO  COMPUTE  THE  BRACKET  LONGITUDE 

21 

c* 

VALUES  IS  A  SIMPLE  ALGEBRAIC  ALGORITHM  INVOKING 

22 

c* 

INTEGER  ALGEBRA.  SEPARATE  EOUATIONS  ARE  REQUIRED 

23 

c* 

DEPENDING  ON  WHETHER  VEST  <POS>  OR  EAST  (NEG) 

24 

c« 

LONGITUDES  ARE  USED. 

25 

c* 

26 

c* 

LIMITATIONS- 

ANY  LONGITUDE  VALUE  FROM  0.00  TO  179.99  CAN  BE 

27 

c* 

USED  AS  INPUT  BUT  THE  BRACKET  LONGITUDE  VALUES 

20 

c* 

WILL  ALWAYS  BE  30.0  TIMES  SOME  INTEGER  NUMBER 

29 

c* 

( Ot  0 i  30 • 0  *  6  0 • 0 •  ETC*)* 

30 

c* 

31 

c* 

IN  THE  CASE  WHERE  THE  INPUT  VALUE  IS  A  BRACKET 

32 

c* 

LONGITUDE  AND  A  POSITIVE  NUMBER.  THE  RETURNED 

33 

c* 

LONGITUDE  BRACKET  VALUES  WILL  BE  THE  INPUT 

34 

c* 

VALUE  PLUS  30.0.  IF  THE  INPUT  VALUE  IS  A  BRACKET 

35 

c* 

VALUE  AND  A  NEGATIVE  NUMBER.  THE  RETURNED  LONGITUDE 

36 

c* 

BRACKET  VALUES  WILL  BE  THE  INPUT  VALUE  AND  THE 

37 

c» 

INPUT  VALUE  MINUS  30.0.  IF  PLUS/MINUS  180.0  IS  THE 

38 

c* 

INPUT  VALUE,  THE  RETURNED  LONGITUDE  BRACKET  VALUES 

39 

c» 

ARE  FICTITIOUS  (PLUS/MINUS  100.0  AND  PLUS/MINUS 

40 

c* 

210.01. 

4  1 

c« 

42 

c* 

INPUT- 

XLNG  =  LONGITUDE  FOR  WHICH  BRACKETING  LONGITUDE 

43 

c* 

VALUES  ARc  DESIRED. 

44 

c* 

45 

c* 

OUTPUT- 

XLNGLO  =  LOWER  BRACKETING  LONGITUDE  VALUE 

46 

Q* 

47 

c* 

XLNGHI  a  UPPER  BRACKETING  LONGITUDE  VALUE 

40 

c* 

49 

c* 

SYSTEM  SUB- 

50 

c* 

PROGRAMS 

51 

c* 

USED- 

I  F  I  X 

52 

c  ♦ 

53 

c* 

USER  SUB- 

54 

c* 

PROGRAMS 

55 

c* 

USED- 

NONE 

56 

c* 

57 

c* 

ESTIMATED 

58 

c  ♦ 

CPU  TIME- 

TO  BE  SUPPLIED 

59 

c* 

60 

c* 

STORAGE 

61 

c* 

REQUIREMENTS 

-  PROCEDURE  PLUS  DATA  OCCUPY  60  WORDS  OF  STORAGE. 

62 

c* 

63 

c* 

program 

64 

c* 

upoates- 

NONE 

65 

c* 

66 

c 

>*****♦*•***•♦* 

67 

c 

68 

DATA  0ELTL0/30.0/ 

r 


'liw11 


69 

C 

ro 

C 

branch  to  positive  or  negative  longitude  computation 

7  1 

C 

72 

IF  <  XLNG )  110.105.105 

73 

c 

74 

c 

LOVER  LONGITUDE  LIMIT  FOR  POSITIVE 

longitude. 

75 

c 

76 

105 

XLNGLO  *  I F I X  ( XLNG/DELTLO >  ♦  DELTLO 

77 

c 

78 

c 

UPPER  LONGITUDE  LIMIT  FOR  POSITIVE 

LONGITUDE. 

79 

c 

80 

XLNGHI  =  XLNGLO  ♦  DELTLO 

81 

GO  TO  120 

82 

c 

83 

c 

LOVER  LONGITUDE  LIMIT  FOR  NEGATIVE 

LONG  I TUDE . 

84 

c 

85 

110 

XLNGLO  =  (  IFIX1XLNG/DELTL0)  -  1>  * 

DELTLO 

86 

c 

87 

c 

UPPER  LONGITUDE  LIMIT  FOR  NEGATIVE 

LONGITUDE. 

68 

c 

89 

XLNGHI  *  XLNGLO  ♦  DELTLO 

90 

c 

91 

c 

TERMINATION. 

92 

c 

93 

120 

RETURN 

94 

END 

58 


■Mini  r--*— 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 
1  1 
12 

13 

14 

15 

16 
17 
lb 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 


CBRLAT  BRACKET  LATITuOE/P.  L «  HER0D/02  JAN  1979 
C 

SUBROUTINE  BRLAT  ( XL AT  >  XLATLO <  XLATHI) 


C 

C*** 
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c 


♦♦♦A**************************************************************** 


program  id- 
met  analyst- 
STS  ANALYST- 
PROGRAMMER- 

CREATED  ON- 

OESCRIPTION- 


mETHOD- 


limitations- 


input- 


OUTPUT- 


SYSTEM  SUB¬ 
PROGRAMS 
USED- 

USER  SUB¬ 
PROGRAMS 
USED- 

ESTIMATED 
CPU  TIME- 


BRLAT 

CAPT  PATRICK 
MAJ  ROGER  C. 
CAPT  PATRICK 

02  JAN  1979 


L.  HER00,  USAFETAC/DnS.  EXT  5412 
WH1T0N,  USAFETAC/DNS,  EXT  5412 
L.  HEROD.  USAFETAC/DNS.  EXT  5412 

PROJECT  192301 


THIS  SUBROUTINE  SUBPROGRAM  FINDS  THE  LATITUDINAL 

values  In  tenths  of  degrees  that  bracket  a  given 

LATITUDE.  THE  LATITUDE  BRACKET  VALUES  ARE  IN 
INCREMENTS  OF  15.0.  I.E.,  A  GIVEN  LATITUDE  OF  27.46 
HAS  BRACKET  VALUES  OF  15.0  AND  30.0. 

THE  METHOD  USED  TO  COMPUTE  THF  BRACKET  LATITUDE 
VALUES  IS  A  SIMPLE  ALGEBRAIC  ALGORITHM  INVOKING 
INTEGER  algebra,  separate  EQUATIONS  ARE  REQUIRED 
DEPENDING  ON  WHETHER  NORTH  ( POS )  OR  SOUTH  ( NEG  > 
LATITUDES  ARE  USED. 

ANY  LATITUDE  VALUE  FROM  0.00  TO  PLUS/MINUS  89.99 
CAN  BE  USED  AS  INPUT  BUT  THE  BRACKET  LATITUDE 
VALUES  WILL  BE  15.0  TIMES  SOME  INTEGER  NUMBER 
<0.0.  15.0.  30.0.  45.0.  ETC.).  IN  THE  CASE  WHERE 
THE  INPUT  VALUE  IS  A  BRACKET  VALUE  AND  A  POS¬ 
ITIVE  NUMBER,  THE  RETURNED  LATITUDE  BRACKET 
VALUES  WILL  BE  THE  INPUT  VALUE  AND  THE  INPUT 
VALUE  PLUS  15.0.  IF  THE  INPUT  VALUE  IS  A  BRACKET 
VALUE  ANO  A  NEGATIVE  NUMBER.  THE  RETURNED  LATITUDE 
BRACKET  VALUES  WILL  BE  THE  INPUT  VALUE  AND  THE 
INPUT  VALUE  MINUS  15.0.  IF  PLUS/MINUS  90.0  IS  THE 
INPUT  VALUE.  THE  RETURNED  LATITUDE  BRACKET 
VALUES  ARE  FICTITIOUS  (PLUS/MINUS  90.0  AND 
PLUS/MINUS  105.0). 

XLAT  =  LATITUDE  FOR  WHICH  BRACKETING  LATITUDE 
VALUES  ARE  DESIRED. 


XLATLO 
XLATHI 

I  F  I  X 

NONE 

TO  BE  SUPPLIED 


LOWER  BRACKETING  LATITUDE  VALUE 
UPPER  BRACKETING  LATITUDE  VALUE 


STORAGE 

REQUIREMENTS-  PROCEDURE  PLUS  DATA  OCCUPY  58  WORDS  OF  STORAGE, 


PROGRAM 

UPDATES- 


NONE 


59 


DATA  DELTLA/ 15. 0/ 


n  n  non  n  n  n  n  n  n  n  n  n  non 


BRANCH  TO  POSITIVE  OR  NEGATIVE  LATITUDE  COMPUTATION. 
IF  ( XL  A  T )  110.10S.105 

LOVER  LATITUDE  LIMIT  FOR  POSITIVE  LATITUDE. 

I0S  XLATLO  x  IFIX  <  XLAT/DELTLA )  *  DELTLA 

UPPER  LATITUDE  LIMIT  FOR  POSITIVE  LATITUDE. 

XLATHI  =  XLATLO  ♦  DELTLA 
GO  TO  120 

UPPER  LATITUDE  LIMIT  FOR  NEGATIVE  LATITUDE. 

110  XLATHI  =  IFIX  (XLAT/DELTLA)  *  oELTLA 

COVER  LATITUDE  LIMIT  FOR  NEGATIVE  LATITUDE. 

XLATLO  =  XLATHI  -  DELTLA 

TERMINATION. 

C 

120  return 

END 
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